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Abstract 

Using the example of the tidal stream of the Milky Way globular cluster Palomar5 (Pal 5), we 
demonstrate how observational data on tidal streams can be efficiently reduced in dimensionality 
and modeled in a Bayesian framework. Our approach combines detection of stream overdensities 
by a Difference-of-Gaussians process with fast streakline models of globular cluster streams and a 
continuous likelihood function built from these models. Inference is performed with Markov chain 
Monte Carlo. By generating « 10 7 model streams, we show that the unique geometry of the Pal 5 
debris yields powerful constraints on the solar position and motion, the Milky Way and Pal 5 itself. 

All 10 model parameters were allowed to vary over large ranges without additional prior information. 

Using only readily-available SDSS data and a few radial velocities from the literature, we find that 
the distance of the Sun from the Galactic Center is 8.30 db 0.25 kpc, and the transverse velocity is 
253 db 16kms -1 . Both estimates are in excellent agreement with independent measurements of these 
two quantities. Assuming a standard disk and bulge model, we determine the Galactic mass within 
Pal5’s apogalactic radius of 19kpc to be (2.1 d= 0.4) x 10 11 M©. Moreover, we find the potential 
of the dark halo with a flattening of q z = 0.95^o!i2 t° be essentially spherical - at least within 
the radial range that is effectively probed by Pal 5. We also determine Pal5’s mass, distance and 
proper motion independently from other methods, which enables us to perform vital cross-checks. Our 
inferred heliocentric distance of Pal 5 is 23.6^o’y kpc, in perfect agreement with, and more precise than 
estimates from isochrone fitting of deep HST imaging data. We conclude that finding and modeling 
more globular cluster streams is an efficient way for mapping out the structure of our Galaxy to high 
precision. With more observational data and by using additional prior information, the precision of 
this mapping can be significantly increased. 

Subject headings: dark matter — Galaxy: halo — Galaxy: structure — Galaxy: fundamental parame¬ 
ters — Galaxy: kinematics and dynamics — globular clusters: individual (Palomar 

5 ) 


1. INTRODUCTION 

Every satellite orbiting a host galaxy produces tidal 
debris. Yet, we find debris only around very few globular 
clusters or dwarf galaxies. This is because the amount 
and shape of tidal debris differs substantially from satel¬ 
lite to satellite. That is, every satellite leaves its own 
unique trace in the gravitational potential of its host. 
This being said, we can expect tidal debris to contain 
substantial information about the satellite, its orbit, 
and its host. Using the example of the most prominent 
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globular cluster stream discovered so far, Palomar 5, 
we are going to demonstrate in the following how this 
information can be extracted through sophisticated 
stream detection, numerical modeling and Bayesian 
inference. 

With the onset of the era of survey science, detec¬ 
tions of tidal debris have become numerous - as 
extended overdensities in surface density; first from 
photographic plates (Grillmair et al. 1995; Kharchenko 
et al. 1997; Lehmann & Scholz 1997; Leon et al. 2000; 
Lee et al. 2003), and more recently in the data from 
the Sloan Digital Sky Survey (Odenkirchen et al. 2001, 
2003; Grillmair 2006; Grillmair & Dionatos 2006b; 
Grillmair & Johnson 2006; Grillmair & Dionatos 2006a; 
Belokurov et al. 2006a, 2007). This number grew as the 
homogeneity and coverage increased with subsequent 
data releases (Grillmair 2009; Balbinot et al. 2011; 
Bonaca et al. 2012; Grillmair et al. 2013; Bernard et al. 
2014; Grillmair 2014). Newly available imaging surveys 
like ATLAS (Koposov et al. 2014) or PAndAS (Martin 
et al. 2014) add to this growing list of tidal features. 
Some tidal streams were also, or solely, found as coherent 
structures in velocity space (Helmi et al. 1999; Newberg 
et al. 2009; Williams et al. 2011; Carlin et al. 2012; 
Sesar et al. 2013; Sohn et al. 2014). With the richness 
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of ongoing and future surveys like APOGEE, DES, 
LAMOST, PanSTARRS, RAVE, SkyMapper, Gaia, and 
LSST (Perryman et al. 2001; Steinmetz et al. 2006; 
Keller et al. 2007; Ivezic et al. 2008; Kaiser et al. 2010; 
Eisenstein et al. 2011; Rossetto et al. 2011; Cui et al. 
2012; Zhao et al. 2012), we will be supplied with ever 
growing amounts of observational data on tidal debris 
of known and unknown satellites within the Milky Way 
disk and in the halo. 

Despite this amount of data - available or yet to 
come - few attempts have been made to exploit the 
wealth of information encoded in these coherent struc¬ 
tures. Most investigations studied theoretical cases 
(e.g. Murali & Dubinski 1999; Eyre & Binney 2009b,a; 
Varghese et al. 2011; Peharrubia et al. 2012; Sanders 
2014), and only a few have actually tried to model the 
observations directly (e.g. Dehnen et al. 2004; Fellhauer 
et al. 2007; Besla et al. 2010; Koposov et al. 2010; Lux 
et al. 2013), as appropriate modeling techniques are 
only becoming available now. The two most notewor¬ 
thy examples for streams that have been modeled by 
utilizing the observational data self-consistently are the 
Sagittarius stream and the GD-1 stream, which will be 
discussed in Sec. 4 in detail. The modeling of these two 
tidal streams have demonstrated the power (but also 
the caveats) of using tidal streams to infer the mass and 
shape of the Milky Way’s gravitational potential. 

The reason for the relatively low number of fully 
modeled tidal streams is twofold: first of all, generating 
a realistic model of a satellite’s debris usually requires 
Wbody simulations of some sort, which are computa¬ 
tionally expensive and therefore do not allow for large 
model-parameter scans. However, a realistic stream 
model is necessary, since it has been shown that streams 
do not strictly follow progenitor orbits (Johnston 1998; 
Helmi et al. 1999), and that fitting a stream with a single 
orbit can lead to significant biases in potential estimates 
(Lux et al. 2013; Sanders & Binney 2013). Angle-action 
approaches for generating stream realizations are a 
promising way of reducing the computational needs 
(Sanders 2014; Bovy 2014). Here we will make use 
of restricted three-body models, so-called streaklines , 
that allow us to generate a realistic stream model at 
small computational cost (Kiipper et al. 2012). The 
second reason is that comparisons of observational data 
to theoretical models are not straightforward. This is 
in most cases due to the faintness of the structures, 
that is, due to the heavy foreground and background 
contaminations from the Galaxy itself. Probabilistic ap¬ 
proaches to such a problem, however, require a manifold 
of stream models, which is why most stream modeling 
approaches so far are only applied to (uncontaminated 
and error-free) simulations of satellite streams (Lux 
et al. 2013; Bonaca et al. 2014; Deg & Widrow 2014) 
or remain qualitative (Belokurov et al. 2006b; Fellhauer 
et al. 2006; Lux et al. 2012; Zhang et al. 2012; Vera-Ciro 
& Helmi 2013). Newly available and computationally 
inexpensive generative models finally allow us to make 
such probabilistic statements on the model parameters 
describing observed tidal streams. 

The problems mentioned above can be alleviated, 


if both observational data and stream models are 
simplified in a sensible way, for example, when indi¬ 
vidual, definite members of a satellite’s debris can be 
identified and when their phase-space coordinates can 
be measured with high precision (e.g., RR Lyraes or 
blue horizontal branch stars; Vivas & Zinn 2006; Abbas 
et al. 2014; Belokurov et al. 2014). With the REWINDER 
method, Price-Whelan & Johnston (2013) demonstrated 
that, in this case, orbits of stream stars can be traced 
back in time to the point where they escaped from the 
progenitor. With only a few stream stars the authors 
demonstrated that they are able to tightly constrain 
their model parameters (Price-Whelan et al. 2014). 
However, such an identification of stream members may 
often not be possible due to the sparsity of tracers or due 
to the unavailability of accurate phase-space information. 

In this work, we present an alternative way to ad¬ 
dress both of these problems, i.e. that allows us to 
quickly generate thousands of stream models, and to 
compare them to well-established parts of an otherwise 
heavily contaminated stream in the halo of the Milky 
Way even if large parts of the phase-space information is 
missing. We will make use of a new method for detecting 
stream over densities, computationally efficient streakline 
models (Kiipper et al. 2012), and the FAST-FORWARD 
method, a Bayesian framework developed in Bonaca 
et al. (2014) for constraining galactic potentials by 
forward modeling of tidal streams (see Sec. 2.3). 

Here we are going to apply a variant of the FAST-FORWARD 
method to the tidal debris of the low-mass globular 
cluster Palomar5 (Pal 5; see Sec. 2.1 for a detailed 
description). Due to its unique position high above the 
Galactic bulge and due to its current orbital phase of 
being very close to its apogalacticon, the tidal tails of 
this cluster are so prominent in SDSS that Pal 5 can 
be considered to be the poster child of all globular 
cluster streams. Its tails stand out so prominently 
from the foreground/background that substructure can 
be seen within the stream with high significance (see 
Sec. 2.1.1, and Carlberg et al. 2012). Using a Difference- 
of-Gaussians process, we will detect overdensities, 
measure their positions and extents, and then use them 
in our modeling to compare our models to these firmly 
established parts of the Pal 5 stream (Sec. 2.6). 

Understanding the origin of such stream over densi¬ 
ties may give us detailed insight to a satellite’s past 
dynamical evolution. Combes et al. (1999) found that 
tidal shocks from the Galactic disk or bulge shocks could 
introduce substructure in dynamically cold tidal tails. 
However, Dehnen et al. (2004) used Wbody simulations 
of Pal 5 to show that disk or bulge shocks are not the 
origin of the observed over densities. Instead, the au¬ 
thors suggested that larger perturbations such as giant 
molecular clouds or spiral arms may have caused this 
substructure. In subsequent publications, several other 
scenarios for substructure formation within tidal streams 
have been discussed. A promising source of frequent and 
strong tidal shocks should be given through dark matter 
subhalos orbiting in the Galactic halo (Ibata et al. 2002; 
Johnston et al. 2002; Siegal-Gaskins & Valluri 2008; 
Carlberg 2009). Encounters with these subhalos should 
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be the more frequent the lower the subhalo mass - 
and the hope is that, with better observational data, 
tidal streams may be used as detectors for these dark 
structures (Yoon et al. 2011; Carlberg 2012; Erkal & 
Belokurov 2014). In a promising first attempt, Carlberg 
et al. (2012) quantified substructures in Pal5’s tidal tails 
and compared the number of detected gaps in the stream 
with expectation of Y-body simulations. The authors 
concluded that some structures in the Pal 5 stream are 
likely to be of such an origin. A different approach was 
taken by Comparetta & Quillen (2011), who suggested 
that Jeans instabilities could form periodic patterns in 
cold tidal streams. However, Schneider & Moore (2011) 
showed that tidal streams are gravitationally stable, 
and hence that Pal5’s overdensities cannot be of this 
origin. As an alternative to tidal shocks or gravitational 
instabilities, Kiipper et al. (2008b) and Just et al. 
(2009) found that epicyclic motion of stars along tidal 
tails can create periodic density variations in the tails 
when they are generated assuming a continuos flow of 
stars through the Lagrange points (see also Sec. 2.3.1). 
Kiipper et al. (2010b) and Kiipper et al. (2012) showed 
that these regular patterns of overdensities and un¬ 
derdensities are also present in tidal tails of clusters 
on eccentric orbits about their host galaxy, and that 
they are most pronounced when the cluster is in the 
apocenter of its orbit. Moreover, Mastrobuono-Battisti 
et al. (2012) demonstrated with Y-body simulations 
that clusters on Pal 5-like orbits can produce regularly 
spaced overdensities in their tails. In the following, 
we are going to show with our modeling of the Pal 5 
stream that several of the observed overdensities closer 
to the cluster are very likely to be due to such epicyclic 
motion. These overdensities help us to find the correct 
model of Pal 5 as they stick out more significantly from 
the foreground/background than the rest of the stream. 

This paper is organized as follows: in the next 
Section we are going to introduce our methodology, 
that is, we will describe the observational data this 
work is based on (Sec. 2.1), and explain how we extract 
the locations of stream overdensities from this data 
(Sec. 2.1.1), as we are going to fit our models to these 
overdensities. The modeling procedure for tidal streams 
is presented in Sec. 2.3, where we describe how we 
generate streakline models of Pal5’s tidal tails, how we 
compare the models to the observational data using a 
Bayesian approach, and how we scan the parameter 
space for the most probable model parameters using 
the Markov chain Monte Carlo method. In Sec. 3, we 
show the results of this modeling, first for an analog 
stream generated from an Y-body simulation (Sec. 3.1), 
and then for the real observations of Pal 5 (Sec. 3.2). 
The last part of this paper is a discussion of the most 
probable parameters resulting from our modeling, with 
a comparison to results from other studies (Sec. 4), 
followed by summary and conclusions in Sec. 5. 

2. METHOD 

2.1. Pal 5 observational data 

Being one of the faintest and most extended globular 
clusters of the Milky Way, Palomar 5 was first discovered 
on photographic plates from the Palomar Observatory 
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Fig. 1.— Probability-density contours of Pal 5 stars from our 
matched-filter analysis of the SDSS data set. Pal 5 lies at (lcos (b), 
b) of (0.59, 45.86) deg and has been masked out (red cross). The 
leading tail is to the right of the cluster, the trailing to the left. 
The large structure at (2.64, 46.80) deg is the residual of the fore¬ 
ground cluster M 5. Upper panel: density contours smoothed with 
a Gaussian kernel with a FWHM of 1.8 deg showing the large scale 
variations across the stream. Lower panel: same as upper panel 
but for a FWHM of 0.23 deg. Smaller kernel sizes reveal more 
substructure but also generate more statistical artifacts. For our 
analysis we subtracted the upper panel from the lower panel (see 
Fig. 3). 

Sky Survey (Abell 1955). Years later and to the great 
benefit of tidal stream research, it happened to lie within 
the footprint of the Sloan Digital Sky Survey (SDSS), 
although unfortunately quite close to the southern edge 
of the survey footprint. Using a matched filter on SDSS 
commissioning data, Odenkirchen et al. (2001) detected 
the tidal tails of Pal 5 stretching 2.6 deg on the sky. This 
marked the discovery of the first coherently mapped 
tidal stream. With subsequent data releases, the tails 
were found to span 22 deg within the SDSS footprint 
(Rockosi et al. 2002; Odenkirchen et al. 2003; Grillmair 
& Dionatos 2006a), where the trailing tail can be traced 
for about 18.5 deg, while the leading tail is cut off by the 
edge of the survey at a length of only 3.5 deg. Successive 
attempts to detect the trailing tail beyond this extent, 
e.g., by using neural networks (Zou et al. 2009), had not 
been successful. However, with SDSS Data Release 8 the 
photometric uniformity and coverage of the Pal 5 stream 
increased (Aihara et al. 2011). In this data, Carlberg 
et al. (2012) found the trailing stream to extend out to 
23.2 deg from the cluster center, where the stream runs 
off the footprint. 

To produce a density map of the Pal 5 stream for 
our analysis, we apply the method outlined in Balbinot 
et al. (2011) to the SDSS DR9 data (Ahn et al. 2012). 
We select objects brighter than g = 22.5 and classified as 
stars by the SDSS pipeline. Additional color cuts were 
performed following the same recipe as in Odenkirchen 
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Fig. 2. — Histogram of significances computed for each pixel in 
our residual map using Eq. 1. The distribution of significances 
is approximately normal-distributed. Overdensities in the resid¬ 
ual map are detected by searching for regions with more than 10 
connected pixels with significance > 3cr. 

et al. (2003). The resulting map is represented on a 
grid of 0.03 x 0.03 deg bins. In Fig. 1 we show this map 
smoothed with a Gaussian kernel, using two different 
kernel widths (0.23 and 1.8deg, respectively). In the 
lower panel, Pal5’s tidal tails can be clearly seen as they 
emanate from the cluster at (lcos(b), b) = (1, 46) deg 
and extend to about (-3, 46) deg for the leading tail, 
and about (15, 39) deg for the trailing tail. Both tails 
are limited by the SDSS footprint, but the trailing tail 
gets significantly fainter at the far end due to a wider 
dispersion of the stream, and also due to the larger 
extinction and contamination from the bulge at lower 
Galactic latitudes. 


As can be seen in the lower panel of Fig. 1, the 
foreground/background around the stream shows many 
density enhancements of Pal 5-like stars. The most 
obvious of these features is the residual of the fore¬ 
ground cluster M5 at (3, 47) deg. In order to distinguish 
between parts that belong to the Pal 5 stream and 
contaminations, we are going to extract regions within 
the SDSS footprint in which Pal 5-like stars are statis¬ 
tically over-dense compared to foreground/background 
fluctuations. 


2.1.1. Overdensity extraction 

Carlberg et al. (2012) detected and characterized 
statistically significant density variations within the 
Pal 5 stream, which had already been noticed by 
Odenkirchen et al. (2003). Here we are going to make 
use of such local density enhancements along the stream 
by exploiting the fact that these regions are the most 
likely places to find Pal 5-like stars, and hence we will 
require our models to have high densities in these regions. 

To enhance a certain scale of density fluctuations, 
which are of the order of the stream width (^ 0.2 deg, 
Carlberg et al. 2012), while at the same time remov¬ 
ing foreground/background variations and large-scale 
variations of the tail density itself, we adopted the 
Difference-of-Gaussians process. This process is a 
simplification of the “Laplacian of the Gaussian” tech¬ 


TABLE 1 
Overdensities 


Name 

lcos(b) [deg] 

b [deg] 

Significance 

FWHM [deg] 

T19 

14.66 

39.19 

3.48 

0.15 

T18 

13.71 

39.64 

3.41 

0.14 

T17 

12.96 

40.06 

3.96 

0.28 

T16 

12.28 

40.48 

4.38 

0.23 

T15 

11.77 

40.84 

4.36 

0.34 

T14 

11.28 

41.25 

4.03 

0.23 

T13 

9.65 

42.33 

4.34 

0.21 

T12 

8.85 

42.68 

3.61 

0.16 

Til 

8.39 

42.75 

3.48 

0.17 

T10 

7.26 

43.69 

4.23 

0.24 

T9 

6.02 

44.37 

3.93 

0.19 

T8 

5.43 

44.68 

3.76 

0.19 

T7 

4.74 

44.89 

5.29 

0.33 

T6 

4.14 

45.11 

7.43 

0.74 

T5 

3.70 

45.15 

5.22 

0.31 

T4 

3.04 

45.45 

10.31 

0.57 

T3 

2.50 

45.54 

4.17 

0.22 

T2 

1.82 

45.80 

5.06 

0.42 

T1 

1.27 

45.83 

6.67 

0.48 

LI 

-0.06 

45.79 

7.94 

0.68 

L2 

-0.95 

45.89 

7.93 

0.64 

L3 

-1.37 

45.82 

3.29 

0.11 

L4 

-1.73 

45.84 

3.93 

0.21 

L5 

-2.05 

46.13 

3.16 

0.11 


nique employed in computer vision for blob detection 
(e.g. Lindeberg 1998). First, we masked out the central 
part of the cluster before convolving the map with the 
Gaussians. Then, we created a “small-feature map” by 
convolving the pixel values of the matched-filter, E, 
with a normal function, A/", of a certain width, or, where 
g i = 0.115 deg in our case (lower panel of Fig. 1). From 
this map we subtracted a “background map” smoothed 
with a larger kernel size, cr 2 ~ 8 o\ (upper panel of 
Fig. 1). To compute the significances for each pixel in 
this “residual map”, we followed Koposov et al. (2008) 
and defined the significance, S', for each pixel to be 


S = V47T0T 


E*(A/>i)-A/> 2 )) 

yjT, *Af((T2) 


(1) 


Fig. 2 shows a histogram of significances, S, of the pixels 
of our residual map. As Koposov et al. (2008) point 
out, the Poisson distribution of sources should yield 
a normal distribution for S with a variance of 1. We 
consider pixels with values of S larger than 3 as sta¬ 
tistically significant (red part of the histogram in Fig. 2). 


On this flat, background subtracted map, we searched 
for extended overdensities with SExtractor (Bertin 
& Arnouts 1996). The code was configured to search 
for groups of at least 10 connected pixels with S > 3. 
To compute the FWHM and barycenters, the analysis 
threshold was lowered to S > 2. Note that we masked 
the regions within « 1.5 deg of the borders of the 
SDSS footprint in order to avoid edge effects that arise 
from the Difference-of-Gaussians process. Finally, we 
excluded overdensities found near the south-east edge of 
the footprint due to its proximity to the Galactic bulge 
(shaded region in Fig. 3). This edge lies at « 42 deg 
from the Galactic centre, which seems close enough to 
introduce significant noise in the matched filter. 


The resulting 24 overdensities (19 in the trailing 
tail and 5 in the leading tail; shown in Fig. 3 and 
separated by a horizontal line) with a significance above 
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Fig. 3.— Residual map of the region around Pal 5 (red cross) as obtained by subtracting the upper panel of Fig. 1 from the lower panel. 
The greyscale shows the local significance of each pixel as given from Eq. 1. Barycenters of the 24 detected overdensities are marked by 
pairs of vertical lines, and their FWHM are illustrated by disks (see also Tab. 1). We excluded the grey shaded region from our analysis, 
where the contaminations from the Galactic bulge, or border effects produce random features. Contours correspond to significances of [2, 
3.5, 5, ...]. 


TABLE 2 

Radial velocities of probable stream stars 


lcos(b) [deg] 

b [deg] 

V R [kms x ] 

A Vr [kms x ] 

6.19 

44.44 

-57.08 

0.29 

4.90 

44.79 

-54.60 

0.36 

4.70 

44.89 

-66.93 

0.44 

3.13 

45.43 

-53.78 

0.38 

3.11 

45.34 

-56.79 

0.32 

3.05 

45.25 

-51.83 

0.23 

2.44 

45.52 

-55.52 

0.25 

2.34 

45.58 

-54.73 

0.37 

0.47 

45.77 

-58.28 

0.48 

-0.17 

45.68 

-57.91 

0.31 

-0.19 

45.75 

-58.33 

0.46 

-0.47 

45.78 

-59.83 

1.16 

-0.48 

45.81 

-59.25 

0.45 

-0.92 

45.92 

-47.00 

0.27 

-1.22 

45.82 

-60.96 

0.24 

-2.22 

45.98 

-60.83 

0.34 

-2.27 

45.87 

-65.02 

0.90 


3 a are listed in Tab. 1, sorted from large to small values 
of lcos (b). We also give their significances and FWHM 
as determined by Sextractor. 

In our streakline modeling, we will use the loca¬ 
tions and widths of the detected blobs as the positions 
and positional uncertainties of the tidal tail centerline, or 
in other words as regions with highest local probability 
of finding Pal 5 stars. We will use the over densities 
to calculate a likelihood of each model (see Sec. 2.6). 
In addition to these high-probability regions we will 
compare our streakline models to the observed positions 
and radial velocities of stars in the tidal tails of Pal 5. 

2.1.2. Radial velocities 

Kinematic information on stars in a tidal stream can 
help significantly to constrain model parameters (see 
e.g. Pearson et al. 2015). The proper motion of Pal 5 
has to be assumed as being largely unknown, since three 
different photographic-plate measurements all yielded 
different results, which all disagreed within their stated 


uncertainties (Schweitzer et al. 1993; Scholz et al. 1998, 
and see discussion in Odenkirchen et al. 2001). However, 
there is accurate information on its radial velocity and 
even on the radial velocities of some stars within its 
tidal tails. 

Odenkirchen et al. (2002) measured radial veloci¬ 
ties of 17 giants within the tidal radius of Pal 5 using 
VLT/UVES. From their analysis of these confirmed 
cluster members, they determined the heliocentric 
radial velocity of Pal 5 with a very high precision to be 
—58.7 d= 0.2 kms -1 . This is much more precise than, 
e.g., the velocity uncertainty of the Solar motion with 
respect to the Galactic center. For this reason, we will 
neglect this uncertainty of 0.2 km s -1 on Pal5’s radial 
velocity in the subsequent analysis. 

In a later study, Odenkirchen et al. (2009) pre¬ 
sented VLT/UVES radial velocities of further 17 giants 
with locations well outside Pal 5 but, in projection, close 
to the dense parts of its tidal tails. The locations and 
velocities of these stars are listed with their individual 
uncertainties in Tab. 2. A horizontal line indicates, 
which stars lie in the trailing part of the tail (above), 
and which are in the leading (below). Similar to the 
overdensities, we will compare our models to these 17 
positions and velocities in order to compute a model’s 
likelihood (see Sec. 2.6). Note though that the locations 
of the overdensities represent the summarized positions 
of several stream stars, whereas the radial velocities 
represent individual stars in phase space that may not 
be representative of large fractions of stream stars. 

Recently, Kuzma et al. (2015) published 39 additional 
radial velocity measurements of red giants along the 
Pal 5 stream. The radial velocity gradient along the Pal 5 
stream of their larger sample of 1.0 zb 0.1 kms -1 deg -1 
matches the one determined by Odenkirchen et al. 
(2009) of 1.0 zb 0.4 kms -1 deg -1 . Such additional data 
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will be valuable for future modeling of Pal 5. 

2.2. Pal 5-analog simulation data 

In order to test our methodology of overdensity 
extraction and streakline modeling, we ran an 7V-body 
model of a Pal 5-like cluster. The parameters for this 
simulation were taken from the results of our first, 
preliminary modeling of the observations, and hence 
this analog 7V-body model of Pal 5 is sufficiently close 
to reality (i.e. our final fitting results) to allow such 
a comparison. For the integration of the model, we 
used the freely available, and GPU-enabled 7V-body 
code Nbody6 (Aarseth 2003; Nitadori & Aarseth 2012), 
which we slightly modified such that it could integrate 
the cluster evolution in the same galactic potential as we 
used for our streakline modeling of the Pal 5 observations 
(see Sec. 2.4 for a detailed description of the galactic 
potential). 

In all our modeling, we kept the disk and bulge potential 
fixed as described in Sec. 2.4, since we mainly focussed 
on the potential of the dark halo. We assumed an NFW 
potential for the halo (Navarro et al. 1997), which we 
allowed to be flattened along the z axis, i.e. perpendicu¬ 
lar to the galactic disk. For the analog model we chose 
a scale mass of MhclIo — 1.475 x 1O 12 M 0 , and a scale 
radius of rn a io = 36.5 kpc. The flattening was chosen to 
be q z = 0.94. The scale mass and scale radius correspond 
to a virial mass of M 200 = 1.60 x 1O 12 M 0 , a virial 
radius of r 2 oo = 217kpc, and a concentration of c = 5.95. 

We set up the initial configuration for the star cluster 
model using the publicly available code McLuster 9 
(Kiipper et al. 2011a,b). For simplicity, we arbitrarily 
chose the initial cluster to be a Plummer sphere of 50 000 
single stars of 0.44 M 0 each, with a half-mass radius 
of 11 pc. This simple model is sophisticated enough to 
reproduce the shape of the tidal tails. We followed the 
evolution of the cluster for 6Gyr, during which it lost 
10 000 M 0 of stars into the tidal tails, resulting in a final 
(i.e. present-day) cluster mass of « 12OOOM 0 . Hence, 
the mass loss rate of this model is about 1.67M 0 Myr -1 . 
The orbit of the cluster was chosen to be similar to our 
best-fit results from the observational study (cf. Tab. 7), 
and yields an apocentric radius of 19.0 kpc, a pericentric 
radius of 7.4 kpc, and hence an orbital eccentricity of 
0.44. 

We then took a snapshot of the resulting distribu¬ 
tion of stars from the location of the Sun (see Sec. 2.5 
for details) and inserted this stellar distribution into the 
SDSS source catalog with an offset of 4 deg from the 
original Pal 5. Each star particle from the simulation 
was randomly assigned a magnitude and color according 
to a stellar mass drawn from a present-day mass func¬ 
tion (Kroupa 2001). These were generated such that 
the simulation particles roughly match Pal5’s stellar 
population by using a PARSEC isochrone (Bressan et al. 
2012) with an age of 11.5 Gyr (Dotter et al. 2011) and a 
metallicity of [Fe/H] = —1.3 (Smith et al. 2002). 

It turned out that, when we used the full range of 

9 https://github.com/ahwkuepper/mcluster.git 


stellar masses down to 0.08 M 0 , the number of stars in 
the tails of our Pal 5-analog simulation bright enough 
to be within the SDSS magnitude limits is lower than 
the number of observed stars in the tails of the real 
Pal 5. Thus, to get a similar number of stars in the 
simulated tidal tails as in the original Pal 5 stream, 
we truncated the present-day mass function at the 
low-mass end of 0.5 M 0 , so that enough stars fell into 
the SDSS magnitude limits. This may be interpreted 
as a first indication that the real Pal 5 is losing stars at 
a significantly higher rate than our analog model, but 
could also originate from a different degree of orbital 
compression of the two streams. That is, if the real Pal 5 
is closer to apogalacticon or on a more eccentric orbit, 
the stream may be more strongly compressed and hence 
the surface density of stars in the tails may be higher. 

On this new data set we applied the same methodology 
as described above, and extracted in total 23 new 
over dense regions for the analog stream. However, when 
we inserted the stream into the SDSS catalog, we had 
to shift it by 4 deg so that it did not overlap with the 
real Pal 5. Due to the angular shift, our analog stream 
was more affected by the Galactic bulge and ran earlier 
off the SDSS footprint. Hence, the analog model did 
not produce statistically significant overdensities at the 
far end of the trailing tail and is therefore 4 deg shorter 
in one direction. The leading tail, on the other hand, 
extends about 2 more degrees than the observed Pal 5. 
Thus, the two streams, observed and simulated, have a 
similar length and a comparable number of overdense 
regions. 

We also generated a set of mock radial velocities 
from the simulation. From the set of 7V-body particles 
we randomly picked 17 particles with similar locations 
on the sky as the 17 giants Odenkirchen et al. (2009) 
observed. As uncertainty for each radial velocity we 
chose a conservative value of 0.5kms -1 (cf. Tab. 2). 
These are the ingredients that went into our streakline 
modeling. In the following we are going to describe the 
details of this modeling process. 

2.3. Streakline modeling 

Ideally, we would like to create a comprehensive set 
of direct 7V-body models of Pal 5 and compare them to 
the available set of observations. However, even for such 
a sparse and low-mass globular cluster like Pal 5, a full 
7V-body computation of the cluster and its forming tidal 
tails, covering several Gyr of cluster evolution, takes a 
few days on a modern GPU workstation with the high- 
end code Nbody6. A parameter-space study involving 
many hundred thousand realizations of Pal 5 in different 
Galactic potentials and on different orbits, like we en¬ 
visage here, is out of question with such an approach. 
Therefore, we will use streakline models to generate real¬ 
istic globular cluster stream realizations within a matter 
of seconds, where we follow the FAST-FORWARD modeling 
approach developed in Bonaca et al. (2014). The concept 
of this method is described in the following. 

2.3.1. Escape of stars from a cluster 

Kiipper et al. (2012) showed that the formation of tidal 
tails of Pal 5-like clusters on typical orbits within a galaxy 
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can be very well understood and approximated by re¬ 
stricted three-body simulations in which massless test 
particles are integrated together with a massive cluster 
particle within a background galaxy potential. Kiipper 
et al. (2012) demonstrated that cluster stars preferen¬ 
tially escape through the Lagrange points of the cluster- 
galaxy system ( Rl ), which can be calculated to lie at a 
distance of 

_( GM c \ 1/3 

rL \n 2 c -d 2 $/dR 2 J ’ 

from the cluster center along the line connecting the 
galactic center and the cluster center (King 1962). In 
this equation, G is the gravitational constant, M c is the 
cluster mass, TL C is the instantaneous angular velocity of 
the cluster on its orbit within the galaxy, and d 2 <&/dR 2 
is the second derivative of the galactic potential, 4>, with 
respect to the galactocentric radius of the cluster, R c . 
Moreover, Kiipper et al. (2012) showed that for low-mass 

satellites as Pal 5 the velocities of escaping stars, V, are 
preferentially such that they match the angular velocity 
of the cluster center with respect to the galactic center. 
That is, we can generate the initial conditions for the 
test particles ( R , V) by using 

= X (Rc T r L ) =F Sr*, (3) 

1L C 

V = h x (y c =Fft c r L ) tSv*, (4) 

with i denoting the x, y, and z components, and the 
upper sign being for the leading tail while the lower one 
is for the trailing tail. The two random offsets, Sr and 
5v, can be added to the initial conditions of the test 
particles to emulate the scatter around these idealized 
escape conditions as seen in full 7V-body simulations 
(Kiipper et al. 2012; Lane et al. 2012; Bonaca et al. 
2014; Amorisco 2014; Gibbons et al. 2014; Fardal et al. 
2014). 

Bonaca et al. (2014) demonstrated that with a re¬ 
alistic choice for the scatter in the escape conditions it 
is possible to recover the underlying model parameters 
of direct 7V-body simulations of a dissolving globular 
cluster in an analytic galaxy potential. Here we chose 
a slightly different approach from the FAST-FORWARD 
modeling of Bonaca et al. (2014) by setting these random 
components to be zero in order to create the kinemati¬ 
cally coldest stream possible. We do this since we are 
interested in reproducing the epicyclic substructure in 
the stream with as few test particles as possible (to 
save computational time), and this substructure will be 
preferentially produced by the coldest escapers (Kiipper 
et al. 2012). Moreover, Kiipper et al. (2008a) and Just 
et al. (2009) showed that most escapers will evaporate 
from a star cluster with just about the right amount of 
energy to escape. This is the family of escapers that 
will produce the substructure in the stream, whereas 
higher-temperature escapers or ejecters (e.g. from 
three-body interactions) will basically only introduce 
noise in the substructure pattern or not orbit within 
the tails at all. In Sec. 3.1.2 we will demonstrate that 


this approach recovers the model parameters of direct 
Ai-body simulations of a dissolving globular cluster on a 
Pal 5-like orbit perfectly well. 

2.3.2. Recapture of stars and the edge radius 

This picture of the escape process, which we use to 
generate our streakline models, is somewhat different 
from the notion introduced by King (1962), who as¬ 
sumed that cluster stars are stripped at perigalacticon 
such that the limiting radii of the globular clusters we 
observe today correspond to their tidal radii (Eq. 2) at 
the perigalactica of their orbits. However, it was shown 
with 7V-body simulations in Kiipper et al. (2010a) and 
Webb et al. (2013) that clusters adjust to the mean tidal 
field of an orbit rather than to the perigalactic tidal field. 
Kiipper et al. (2012) demonstrated that the tidal radius 
according to King’s formula for a cluster on an eccentric 
orbit gets so small at perigalacticon that stars cannot 
escape from there 10 . Instead, they will get recaptured 
by the star cluster while it is moving into apogalacticon, 
as the tidal radius quickly grows to its maximum value. 
This demonstrates the importance of the cluster’s 
gravitational acceleration on the trajectories of escaping 
stars. In fact, Kiipper et al. (2012) showed that, even for 
clusters on circular orbits, the trajectories of tail stars 
are significantly affected by the cluster mass. This gets 
more important for more eccentric orbits, since tail stars 
will have periodic strong encounters with the cluster 
in apocenter. For very eccentric orbits, this effect even 
influences the shape of the tidal tails at large distances 
from the cluster, as the compression of the tidal tails in 
apocenter is so significant that distant debris is pushed 
back into the cluster vicinity and remains close to the 
cluster for a significant fraction of the orbital period. 

Kiipper et al. (2012) suggest to introduce a lower 
limit on the tidal radius from Eq. 2 in order to prevent 
recapture of test particles. In our investigation we chose 
this edge radius to be at least equal to Pal5’s half-light 
radius of 20 pc (Harris 1996), and iteratively increase it 
for a specific test particle if it gets recaptured by the 
cluster. The cluster itself is represented by a softened 
point mass, M c , with a softening length of the order of 
the cluster’s half-light radius of 20 pc. The influence 
of this latter number on the actual shape of the tails is 
negligible, which is why we fixed it to some reasonable 
value instead of making it a free parameter in the 
modeling process. For more massive and much more 
extended systems like the Sagittarius dwarf galaxy, the 
progenitor should be modeled with more care as the 
extent of the system more strongly influences the shape 
of the tidal tails (Gibbons et al. 2014). 

2.3.3. Test-particle integration 

The streakline models are generated in an analytic 
background galaxy potential, <f>, which will be specified 
in the next Section. Given the more or less well 
constrained 6D phase-space coordinates of Pal 5 at 
the present day, we integrate the cluster particle from 
this initial position backwards in time for a specified 
integration time, using a 4th-order Runge-Kutta 

10 Unless they have high excess energies, e.g. from three-body 
encounters or very strong tidal shocks. 



integrator. From this point in the past (—U n t) we 
integrate back to the future until we reach the present 
day, while releasing test particles at a fixed time interval, 
dt re i ea se • Once a test particle is released from the cluster 
particle at a given point in time in the past, the test 
particle is integrated together with the cluster particle 
within the galaxy potential until the present day. This 
is repeated K = t int / dt re i ease times until the cluster 
particle has reached its present-day position. The 
test particles do not interact with each other, as the 
inter-particle distance within the stream is too large to 
cause relevant interactions. Thus, each streakline model 
is a combination of K restricted three-body integrations, 
where the longest integrations are over a time interval 
ti nt and the shortest ones are only of length dt re i ease . 
We chose a interval length of dt re i ease = lOMyr, which 
we found to produce well-resolved stream models with 
acceptable computational effort. Pal 5 loses stars at 
a much higher rate, though (probably of the order 
of 20 stars per Myr, when we use our fitting result 
of 7.9M 0 Gyr _1 and assume a mean stellar mass of 
0.4 M 0 ). Our test particles are therefore only tracers 
of the stream, not exact representations of the actual 
stellar densities. 

Note that this approach produces streaklines with 
a constant release rate of test particles. Tidal shocks, 
e.g. when the cluster crosses the galactic disk, do not 
increase the release rate temporarily. This appears to 
be justified as no obvious density variations resulting 
from disk shocks can be found in 7V-body simulations 
of star clusters on Pal 5-like orbits (Dehnen et al. 
2004; Kiipper et al. 2010b). For real clusters and full 
7V-body simulations of star clusters, the eccentricity 
of an orbit does change the mass-loss rate, though, 
as higher eccentricities will make the cluster plunge 
deeper into the central parts of the galactic potential 
(Baumgardt & Makino 2003). However, this increased 
mass loss appears to be spread out over the entire orbit 
of the cluster (Kiipper et al. 2010b). This delayed 
escape of stars that got unbound during a tidal shock 
is due to the relatively long timescale of escape for just 
slightly unbound stars (Fukushige & Heggie 2000; Zotos 
2015), as well as due to the recapture of many stars 
that got temporarily unbound during a tidal shock as 
mentioned above (Sec. 2.3.2). Of course, stars that 
gained a lot of energy during a tidal shock will quickly 
escape from the cluster. This can be seen for example 
in simulations of dwarf galaxies falling radially into the 
potential of a host galaxy. In such extreme cases the 
tidal shocks will produce substructure in the satellite’s 
stream (Penarrubia et al. 2009). However, such stars 
with high excess energies will not contribute to the 
coldest part of the stream that we are interested in. Our 
streakline method therefore produces cold tidal streams, 
which trace the density variations along the streams 
due to orbital compression and stretching of the stream 
in apogalacticon and perigalacticon, respectively, as 
well as due to epicyclic motion. But the density of test 
particles does not directly correspond to the density 
of stars in the stream, as the fraction and distribution 
of kinematically warmer stream particles is not well 
known. The test particles are therefore just tracers of 
the locations where the coldest stream stars pile up. 


Similar approaches of test-particle integrations for 
tidal streams have been used in previous studies 
(Odenkirchen et al. 2003; Varghese et al. 2011; Gibbons 
et al. 2014), but in our case the method is tailored to 
globular clusters and calibrated with 7V-body simula¬ 
tions, which makes the resulting stream models more 
realistic and - as we will show - estimates on, e.g., 
parameters of the Galactic potential more reliable. 

2.4. Galactic potential 

The gravitational potential of the Milky Way (MW) is 
one of the main quantities of interest here. The choice of 
representation or, in most cases, analytical parametriza- 
tion of the potential is not straightforward, and ideally 
it should be kept as flexible as possible. Great flexibility, 
however, implies many degrees of freedom. Given the 
limited amount of data we currently have for Pal 5, we 
want to keep the free parameters at a minimum but still 
be flexible enough to accurately describe the MW po¬ 
tential and get a meaningful comparison to other studies. 


In past studies, the Galactic potential has often 
been parametrized as a three-component model consist¬ 
ing of a bulge, a disk and a halo (e.g. Johnston et al. 
1999; Helmi 2004a; Johnston et al. 2005; Fellhauer et al. 
2007; Koposov et al. 2010; Law & Majewski 2010). We 
will adopt a similar parametrization and focus on the 
halo component, as this is the least constrained, but also 
the most massive component of the Milky Way. The 
other two components, bulge and disk, are kept fixed 
and are modeled as follows. For the bulge potential we 
use a Hernquist spheroid (Hernquist 1990), 


^Bulge (0 


G MBulge 

r + a 


( 5 ) 


where G is the gravitational constant, r = y/x 2 + y 2 + z 2 
is the Galactocentric radius, MBulge = 3.4 x 1O 1O M 0 , 
and a = 0.7 kpc. As disk component we use an axis- 
symmetric Miyamoto & Nagai (1975) potential, 


®Disk{R, z) 


_ GM Disk _ 

\J R 2 + {b + y/z 2 + <?) 2 


( 6 ) 


with R = y/ x 2 + y 2 , Mrnsk = lO n M 0 , b = 6.5kpc, 
and c = 0.26 kpc. These parametrizations and also 
these parameter values for the two stellar components 
are not the most sophisticated or up-to-date choices 
(see e.g. Dehnen & Binney 1998; Binney & Tremaine 
2008; Bovy & Rix 2013; Irrgang et al. 2013; Licquia & 
Newman 2014), but we chose them here for compara¬ 
bility to past studies. In follow-up studies with more 
observational data for Pal 5, and when other constraints 
on the Galactic potential are taken into account, the 
baryonic components should be included in the modeling 
process. However, we argue that the differences between 
the various parametrizations of these components at the 
average Galactocentric distance of Pal 5 are not very 
important, as the halo is the dominant mass component 
at these radii. Moreover, Vesperini & Heggie (1997) 
showed that the Galactic disk has very little effect on 
the evolution of globular clusters at Galactocentric radii 
beyond about 8 kpc, which is where Pal 5 crosses the 
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disk most of the time. 


with results from data sets independent of Pal 5. 


We will use these analytically simple and compu¬ 
tationally inexpensive prescriptions for the bulge and 
disk, and focus on varying the mass, shape and extent of 
the dark halo. As parametrization for the halo we chose 
the NFW density profile (Navarro et al. 1997), as it has 
been shown to be a very good representation for dark 
matter halos in ACDM simulations of structure forma¬ 
tion (see e.g. Bonaca et al. 2014). The gravitational 
potential resulting from this density profile has the form 


$ Halo(r ) = - 


G Mualo 
r 


In 



—)■ 
rHalo J 


( 7 ) 


where MhclIo and rpaio are a scale mass and a scale ra¬ 
dius, respectively, which are left as free parameters in our 
modeling. The scale mass and scale radius can be easily 
related to the more frequently used virial radius, r 2 oo, 
the virial mass, M 2 oo, and the concentration, c. r 2 oo is 
the radius of a sphere with a mean interior density of 
200 times the critical density, p cr it, which has a value of 
Pcrit = 3 H 2 / (8tt G) = 1.26 x 10 -7 M 0 pc -3 when using a 
Hubble constant of H = 67.3kms -1 Mpc -1 (Planck Col¬ 
laboration et al. 2014). The virial mass, that is the mass 
enclosed in a spherical volume of radius r 2 oo, is related 
to the scale mass via 


^200 _ f ^Halo +\ n ( THal ° + r2QQ ^ _ i qA 

Mnalo X^Halo+r 200 \ T Halo J / 

(8) 

and the concentration of the NFW halo is usually 
defined as c = r 2 oo/r Ha io- 


For computational convenience we will use Eq. 7 
for our modeling, but report our results also in terms of 
r 2 oo, Af 2 oo and c to allow comparability to other studies. 
However, due to the flattening along the z-axis, and 
due to the vast extrapolation from Pal5’s apogalactic 
radius of about 19kpc out to a virial radius of about 
200-300 kpc, which solely relies on the assumption of an 
NFW profile for the halo, comparisons between these 
numbers from different studies have to be taken with 
care. 


This is especially important since we also allow 
the potential to be flattened along the z-axis perpen¬ 
dicular to the Galactic disk by introducing a scale 
parameter q z . That is, in the above equation for the 
halo potential we replace r by R? + z 2 /q 2 . Such a 
potential flattening makes the definition of an enclosed 
mass ill-defined. A more reasonable quantity, which 
our modeling may provide, is the acceleration Pal 5 
experiences at its present-day position high above the 
Galactic disk, ap a i 5. This acceleration can be converted 
into a local “circular velocity” 

Vc {rpal 5 ) = \Z\dPal 5 | ^Pal 5, (9) 

and as such allows a convenient comparison to the gravi¬ 
tational field within the Galactic disk, for which accurate 
rotation curves have been measured (e.g. Sofue 2013). 
We will therefore calculate these values from our model¬ 
ing results, and also compute the circular velocity at the 
Solar Circle to compare the performance of our modeling 


2.5. Solar parameters 

One key ingredient for the comparison of model points 
from numerical simulations to data points from obser¬ 
vations is the location and motion of the Sun within 
the Galaxy. We chose our right-handed Galactocentric 
coordinate system such that the Sun lies in the xy-plane 
on the x-axis with the x-axis pointing to the Galactic 
center, and such that its tangential motion is along the 
y-axis in the direction of positive y. 

Observations of the, so-called, S-stars around the 
Galactic center suggest that the distance of the Sun 
from the Galactic center is Rq = 8.33 ± 0.35 kpc 
(Gillessen et al. 2009), which puts it in our coordinate 
system to R Q = (—8.33,0,0) kpc. Using parallaxes 
and proper motions of over 100 masers in the Milky 
Way, Reid et al. (2014) found a similar but somewhat 
more precise value of 8.34 db 0.16 kpc. With this data 
the authors were also able to constrain the tangential 
motion of the Sun within the Galaxy. This motion is 
usually split into the circular velocity of the Galaxy 
at the Solar radius, Vc, and the Solar reflex motion 
(U 0 , U 0 , H 7 ©), that is, the motion of the Sun with re¬ 
spect to an object on a closed circular orbit at Rq. Reid 
et al. (2014) find the total tangential velocity of the Sun 
in y-direction to be V tan = Vc + Vq = 255.2 ±5.1 kms -1 . 

Schonrich (2012) used Hipparcos data to derive 
the Sun’s reflex motion to be 

V Q = (11.1 ± 0.7,12.2 ± 0.5,7.3 ± 0.4) kins -1 , (10) 

which gives us a value of 243kms _1 for the circu- 
lar velocity at the Solar radius. A similar result of 
Vc = 239 zb 5kms -1 was found independently by 
McMillan (2011) using a wider range of observational 
data. 

However, from APOGEE data, Bovy et al. (2012) 
determined a slightly different set of Solar parameters, 
putting the Sun at a distance of 8.1+o.ikpc from the 
Galactic Center with a total tangential velocity of 
Vtan = 242^3° km s _1 , only marginally in agreement 
with the Reid et al. (2014) values. Moreover, while their 
Uq and Wq components are in good agreement with the 
Schonrich (2012) values, they determined the circular 
velocity at the Solar Circle to be as low as 218±6 kms -1 . 

Due to this existing discrepancy, we will not use 
literature constraints on the circular velocity at the 
location of the Sun within the Milky Way nor on its 
tangential velocity. For the coordinate transformations 
between observer coordinates and the Cartesian coordi¬ 
nates of the models we will adopt the Schonrich (2012) 
values for Uq and B 7 ©, but we will keep Vtan, as a free 
parameter. Moreover, we will leave the distance of the 
Sun from the Galactic Center as a free parameter. 

2.6. Comparison of models to observations 

The most important and also most sensitive part for 
drawing conclusions from the modeling of tidal streams 
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is the comparison of the models to the observational 
data. In the following we will use the formalism outlined 
in Bonaca et al. (2014), which is based on the framework 
developed in Hogg et al. (2010) and Hogg (2012). 

Bonaca et al. (2014) defined a likelihood, P, for a 
given set of model parameters, (9, and priors, /, to 
compare stream particles from the Cauda simulation 
(a re-simulation of the Via Lactea II simulation by 
Diemand et al. 2008 including globular cluster streams) 
to particles from streakline models. Similar to this work, 
we will refer to our streakline particles as model points, 
while our data points will be the Pal 5 observations 
or the Pal 5-analog observations, respectively, denoted 
as Xyi • 

Assume we have K model points and N data points. 
We can write the probability, p, of a model point, £&, 
generating a data point, X n , as 

p(X n \x k ,0,I) =J\f(X n \x k ,al + T, 2 k ) , (11) 

where 6 is the set of model parameters that generated 
the model point, and / is any prior knowledge we have 
on the parameters. The model parameters 6 and priors 
/ will be described in detail in Sec. 2.7. We assess 
the probability of each model point to generate a data 
point by either using a 2-dimensional or 3-dimensional 
Gaussian distribution, J\f (X\m,V), for X with mean, 
m, and variance tensor, V, for the overdensities or for 
the radial velocity stars, respectively. That is, for the 
overdensities, the two dimensions are Galactic longitude 
and Galactic latitude, whereas for the radial velocity 
stars we add the velocity as a third dimension. 

The variance tensor V consists of the observational 
uncertainty tensor, cr^, and a smoothing tensor, 

For the overdensities we will use their FWHM as 
positional uncertainty, whereas for the radial velocity 
stars we use a positional uncertainty of 0.115 deg, which 
accounts for the finite width of the stream. Moreover, 
we use their individual uncertainties in velocity from the 
observations as given in Odenkirchen et al. (2009) and 
listed in Tab. 2. The tensor allows us to smooth the 
discrete distribution of points such that their Gaussians 
overlap in phase space. This makes the generative model 
of the stream continuous in phase space, i.e. it turns the 
discrete streakline particles into a probability density 
function in phase space with finite support everywhere. 
The smoothing length in each dimension is kept as a free 
hyper-parameter in 0, over which we will marginalize 
in the end before we report the results for the other 
parameters. 

For each data point, X n , we can determine the 
likelihood of being produced by the whole streakline 
model by marginalizing over the K model points 

K 

p(X n \0,I)=Y / PkP(Xn\k,e,I), (12) 

k =1 

with the normalization factor being P k = 1/K. The 
number of model points K is in principle an arbitrary 
choice, and hence the likelihood should not depend on 


TABLE 3 

Model parameters and range of flat priors 


Parameter 

Prior range 

Description 

MhclIo 

[10 y ; 10 la ] M 0 

Halo scale mass 

r Halo 

[0.1; 100] kpc 

Halo scale radius 

Qz 

[0.2; 1.8] 

Halo flattening 

M Pal5 

[0; 60000J M 0 

Cluster present-day mass 

dM/dt 

[0; 100] M©yr -1 

Average cluster mass-loss rate 

dpal 5 

[20; 27] kpc 

Cluster distance from Sun 

Hu cos((5) 

[—3; —1] masyr -1 

Proper motion in RA 

ns 

[—3; —1] masyr -1 

Proper motion in Dec 

Van 

[212; 292] kms -1 

Solar transverse velocity 

Rq 

[7.5; 9.0] kpc 

Solar Galactocentric distance 


TABLE 4 

Polynomial coefficients for stream centerline 


Coefficient 

Simulation 

Observations 

CiO 

45.9373 

45.9816 

CLl 

-6.13919 x 10 -2 

-2.50988 x 10 -2 

a 2 

-1.3957 x 10 -2 

-2.44554 x 10 -2 

«3 

-1.96972 x 10 -4 

3.32027 x 10 -4 


this number. A model with higher number of particles 
will sample the density distribution of the tidal tails 
better and will reduce noise. However, producing the 
model points is the computationally expensive part of 
the whole modeling process, so it should be kept at a 
reasonable number. Bonaca et al. (2014) find that for a 
computationally efficient scan of the multi-dimensional 
parameter space of the streakline models, we should use 
K being a few times larger than N. While Bonaca et al. 
(2014) use a fixed factor of 5 difference, we will use a 
variable number of model points and different setups for 
the number of data points. However, in any case we will 
keep K/N > 3, and in most cases it is between 5 and 20. 
This variability in K comes from our choice of keeping 
the interval, at which model points are created, fixed to 
dt re i ease = 10 Myr, while leaving the integration time, 
ti n t , as a free hyper-parameter (see next Section). 

To compute the combined likelihood of the set of 
model points creating the set of data points, we assume 
that all data points are independent of each other, so 
that we can multiply their individual likelihoods to yield 

N 

P({X n }\9,I) = ~[[p(X n \9,I). (13) 

n= 1 

This likelihood will be used as input for the MCMC algo¬ 
rithm described in the next Section. Similar applications 
of this formalism can be found in Bonaca et al. (2014) 
and Price-Whelan et al. (2014). 

2.7. Markov chain Monte Carlo (MCMC) machinery 

Using Eq. 13 we can assess the likelihood of a given 
streakline model to produce the observational data. This 
likelihood we can feed into the Markov chain Monte 
Carlo code emcee (Foreman-Mackey et al. 2013) to find 
the posterior probability density for the model param¬ 
eters. emcee is an affine-invariant ensemble sampler, 
which is ideal for massively parallel runs on high- 
performance computing clusters using message-passing 
interfaces such as OpenMP and MPI (e.g., Allison & 
Dunkley 2014). Moreover, it can make use of parallel 
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tempering, which allows us to scan the parameter space 
at different “temperatures” (e.g., Varghese et al. 2011). 
That is, the chains will basically walk through parameter 
space with different stride lengths. Parallel tempering 
is useful for multi-modal likelihood surfaces, as it helps 
walkers to move out of local maxima, and hence enables 
a more efficient scan of the parameter space. 

In total we have 10 free model parameters and 4 
hyper-parameters. The model parameters can be 
grouped into three categories: 

1. halo parameters (scale mass, Mp a io , scale radius, 
rHaio, flattening, q z ), 

2. cluster parameters (present-day mass, Mp a i 5, av¬ 
erage mass-loss rate, dM/dt , distance, dp a i 5, 2 
proper motion components, yU aC os(<5) and Us), 

3. Solar parameters (transverse velocity with respect 
to the Galactic Center, V ta m and Galactocentric 
radius, Rq). 

We set up the model parameters in a very wide range 
around possibly allowed values, which we roughly 
constrain with values from the literature. We will use 
these same bounds for the initial values also as (flat) 
priors, that is, if a walker tries to “walk” outside this 
bound, the likelihood is set to — 00 . The parameters and 
the range of initial values/the priors are listed in Tab. 3. 

Our four hyper-parameters are: the smoothing in 
position for the overdensities, the smoothing in position 
for the radial velocity stars, and the smoothing in 
velocity for the radial velocity stars. The fourth hyper¬ 
parameter is the integration time for the streakline 
models. The latter has to be kept as a free nuisance 
parameter since we do not know a priori how old the 
part of the Pal 5 stream is that we see within the SDSS 
footprint. It may well be that the stream extends 
much further in both leading and trailing direction, or 
that large parts of the tails have dispersed over time. 
Hence, we allow integration times between 2 and 10 
Gyr, and smoothings between 0 and 100 deg or kms -1 , 
respectively. 

We ran emcee on 128 cores of the Yeti cluster at 
Columbia University using a setup of two temperatures, 
with 512 walkers in each temperature. After a “burn-in” 
phase of 200 steps per walker, we checked that the chains 
had converged sufficiently and started the sampling 
phase. For the final posterior probability densities, 
which we are going to report in the following, we ran 
the chains long enough such that each walker had made 
between 500 and 1000 steps. In the following, we are 
going to report the results from the lowest temperature 
chains. 

Since the overdensities and the radial velocity stars are 
different type of data, and since it is not straightforward 
how to combine them into one statistic, we ran 4 different 
analysis approaches on both the analog simulation and 
on the Pal 5 observations. The 4 different configurations 
are the following: 


1. Overdensities : we assess the likelihood of the model 
producing the observational data by only using 
the locations of the overdensities as described in 
Sec. 2.1.1 and listed in Tab. 1 for the observations. 
All over densities are treated the same, and their 
FWHM are treated as the uncertainties of their 
barycenter positions. 

2. Overdensities + radial velocities : in addition to the 
locations of the over densities, we add the locations 
of the radial velocity stars in galactic longitude, 
galactic latitude and radial velocity. 

3. Weighted overdensities : using the FWHM as po¬ 
sitional uncertainties of the overdensities is most 
probably a considerate overestimate. In an attempt 
to assign more weight to the more pronounced (and 
more extended) overdensities, we use their signifi¬ 
cances, 5, as weight factors. The likelihood func¬ 
tion (Eq. 13) then reads 

N 

P({X n }\d,I) = ~[[p(X n \9,I) s . (14) 

n= 1 

We assign the radial velocity stars a significance 
of 1, whereas the overdensities have significances 
ranging from 3 to 10 (Tab. 1). This weighting of the 
overdensities compared to the radial velocity stars 
is somewhat arbitrary, but as it turns out the ac¬ 
tual choice does not have a significant effect on the 
results. The sum of all significances of the 24 over¬ 
densities from the Pal 5 observations is 117 (136 
for the 23 over densities used for the analog simu¬ 
lation), which means that we effectively upweight 
our SDSS data by a factor of 5 compared to the 
radial velocity data. However, this up-weighting is 
not linearly distributed over the over densities, but 
according to their observed significances. Thus, an 
overdensity that is detected, e.g., with a signifi¬ 
cance of 3 and a FWHM of 0.1 deg counts as if 
we detected three Pal 5 stars at the location of the 
overdensity’s barycenter with each having an un¬ 
certainty of 0.1 deg. 

4. Interpolated centerline : as an alternative approach, 
we interpret the locations of the overdensities as 
the centerline (or centroid) of the stream, and in¬ 
terpolate points between these observed overden¬ 
sities. Given Pal5’s convenient shape in Galactic 
longitude, /, and latitude, 5, we fit a third-order 
polynomial of form 

b(l) = ao + a\ l + a2 / 2 + U 3 / 3 (15) 

to the data points using a Marquardt-Levenberg 
algorithm. The resulting coefficients for the ana¬ 
log simulation and for the observations are listed 
in Tab. 4. With this polynomial we generate 117 
points for the observations and 136 points for the 
simulation along the whole extent of the stream. 
We chose these points to be equally spaced in 
Galactic longitude with a spacing width of 0.18 deg 
(0.13deg). The positional uncertainties of these 
centerline points we set to be twice the spacing 
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width, such that we have a continuous representa¬ 
tion of the stream along this fitted centerline. We 
excluded the parts of the centerline between the in¬ 
nermost over densities, i.e. T1 and LI (see Tab. 1). 
In this way, we can use our Bayesian formalism of 
Sec. 2.6 and directly compare the results of this 
centerline fitting to the other approaches. 

Each of the four analysis approaches was run on both 
datasets, the analog simulation and the real observations 
of Pal 5. As mentioned above, we used an MCMC algo¬ 
rithm to explore the model parameter space and deter¬ 
mine the posterior probability densities of each parame¬ 
ter. These densities were determined from the steps of 
the walkers in the Markov chains. We used 512 walk¬ 
ers, and each of them made at least 500 steps, resulting 
in >256000 steps through parameter space. These steps 
provide us with individual histograms for each model pa¬ 
rameter, showing us how often the chains have used a 
certain value for a given parameter. We treat as our 
best-fit parameter the median value of the posterior dis¬ 
tribution for a parameter, and use the values containing 
68% of choices above and below the median as the upper 
and lower uncertainty range, respectively. 

3. RESULTS 

In the next sections, we are going to compare the four 
different analysis approaches described above. First we 
will do this for the analog simulation, which will give us a 
measure of how well the methods recover the simulation 
parameters, i.e. accuracy and precision, and then we will 
report the results for the Pal 5 observations. 

3.1. Pal 5-analog simulation 

As described in Sec. 2.2, we ran a direct 7V-body 
simulation of a star cluster on a Pal 5-like orbit within 
a Milky-Way like potential to test our methodology and 
assess its accuracy. A summary of this assessment is 
shown in Fig. 4 and listed in Tab. 5. The accuracy in 
this plot is defined as (median MCMC value)/(actual 
simulation value). The error bars represent the precision 
of the measurements, that is, the 68% confidence 
intervals. We note that more than 2/3 of the values are 
recovered within these 68% confidence intervals, which 
indicates that we may in fact be overestimating our 
uncertainties. 

All four analysis approaches recover the 10 true 
simulation parameters (and also the 5 derived quantities 
shown in Fig. 4) within their respective uncertainties - 
with only few exceptions, which we will discuss in detail 
in the following. This is very reassuring and shows that 
our methodology works. 

Both the precision and the accuracy of recovering 
the true simulation values vary significantly between 
the model parameters as well as between the methods. 
In general, we find that using only the overdensities 
produces the largest uncertainties (lowest precision). 
This is not surprising since more data also means more 
constraining power. On average, the uncertainties get 
significantly smaller by adding the 17 radial velocity 
stars to our modeling. The exact gain in precision varies 
from parameter to parameter, but lies mostly between 


10-30%. 

The precision can be further improved by adding 
more information, that is, by adding weights to the 
overdensities, or by interpolating between them. The 
run with weighted overdensities and the run with the 
interpolated centerline turn out to achieve similar results 
in terms of precision. The gain in precision compared 
to the run with the overdensities and radial velocities 
can be as large as 50%. From Tab. 5 it is apparent that 
the weighted overdensity analysis approach achieves a 
significantly better accuracy in recovering the simulation 
parameters than the interpolated centerline analysis 
approach. The latter can even introduce such strong 
biases that the true model value lies outside the 68% 
confidence interval. 

In Fig. 5 we show the best-fit model from the weighted 
overdensity approach. The streakline model perfectly 
reproduces the shape of the tails and the radial ve¬ 
locity gradient. The 7V-body simulation created a 
regular overdensity pattern, which originated from 
epicyclic motion within the tails. This pattern has been 
successfully recovered by our Difference-of- Gaussians 
overdensity finder after we inserted the 7V-body model 
into the SDSS field. After the MCMC modeling of 
the detected overdensities, the best-fit streakline model 
nicely reproduces this regular pattern in the vicinity 
of the cluster, and the stream centerline follows the 
detected overdensities further away from the clus¬ 
ter, where the epicyclic over densities are expected to 
be washed out. This demonstrates that epicyclic over¬ 
densities can be used to constrain models of tidal streams. 

The results from the different analysis approaches 
will be discussed in detail, and separately for the 
different model-parameter groups, in the following. 

3.1.1. Halo parameters results (analog simulation) 

All four methods recover the three halo parameters, 
i.e. Mnaio , t Halo and q z , of our analog simulation within 
the uncertainties. The left side of Fig. 4 shows the results 
for the halo parameters for the four analysis approaches. 
As can also be seen from Tab. 5, all parameters are 
recovered within about 20% in terms of accuracy by all 
methods. The best accuracy is achieved by the weighted 
overdensity analysis approach. The interpolated 
centerline method yields a comparable precision, but in¬ 
troduces biases which cause, e.g., the true halo flattening 
parameter to lie just outside the 68% confidence interval. 

The largest uncertainties (lowest precision) we get 
in the scale length of the halo when we only use the 
overdensities, with 1 a uncertainties of +59% and —37%. 
This is expected since we chose the scale length of our 
NFW halo to be 36.5 kpc, which is well outside the 
apogalactic radius of our Pal 5-analog’s orbit. Moreover, 
due to the parametrization of this halo, there is a 
significant degeneracy between the scale length and 
the scale mass of the halo. This can be seen in Fig. 6, 
where we show the posterior parameter distributions 
for the halo parameters for the analysis approach using 
the weighted overdensities. The three blue-shaded 
histograms give the marginalized posterior distributions 
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Fig. 4. — Summary of the accuracies in reproducing parameters from the analog simulation for the four different modeling approaches: 
using only the overdensities, using the overdensities plus the radial velocities, using the weighted overdensities plus the velocities, or using 
an interpolated line going through the overdensities plus the velocities (see Sec. 3 for a detailed description). The accuracy is defined as 
(median MCMC value)/(actual simulation value). The first 10 parameters are the actual modeling parameters for: the halo (Sec. 3.1.1), 
the cluster (Sec. 3.1.2), and the Sun (Sec. 3.1.3), whereas the last 5 parameters are derived quantities (see Tab. 5). 




Fig. 5. — Streakline model (blue semi-transparent points) of the Pal 5-analog simulation using the best-fit parameters from the weighted- 
overdensity analysis approach. The present-day position of the cluster is marked by a red star, whereas as the locations and uncertainties 
of overdensities (upper panel) and radial velocity stars (lower panel) are shown by black markers. 
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TABLE 5 

Accuracy from fitting to the simulation 


Parameter 

Overdensities 

Overdensities 
+ radial velocities 

Weighted 

overdensities 

Interpolated 

centerline 

MhclIo 

r Halo 

qz 

1-19A43 

0.85l°;®? 

1 O 1 !"*” 0 - 25 
i - UO -0.28 

1 1 7+9.50 
± - ±/ -0.46 

1 16+ 0 - 35 
1 ‘ 1O -0.30 

1 10 +0 ' 27 
± - ±U -0.28 

1 02 +(J ' 31 
- L-U -0.28 

1 io +0 - 25 
± - ±U - 0.20 

1 H+0-16 
± - ±i -0.16 

1.08+1^ 

1 1 7+0-23 
± - ± ‘ -0.19 

1 1 Q+0-16 

- Li °-0.17 

Mpal 5 

dM/dt 

dpal 5 
[let COS(A) 

1 1 7+I.23 

^'-o.ss 

9 09+3.53 

z - oz -1.85 

1 00 +0 ‘° 5 
l.UU_o 05 

1 03+ 0 - 12 

1 03+°- 10 

i - u<4 -o.n 

1 1 t^+9-83 
i - iO -0.72 

9 CM+ 4 - 30 
Z.04_2 19 

1 00 +0 ‘° 5 
l.UU_o 05 

o.99l° 0 ;° 0 ? 

0.98l°;°? 

1 19+0-50 

l-iy_ 0 .57 

1 03+ 4 - 75 
l.UO_o 70 

1 nn+°- 04 
l.UU_o 04 

1 nn+°- 07 

l.UU_o 07 

1 nn+ 0 - 07 
4 - uu -0.07 

0 41+0-50 

U -^ 4 -0.24 
n 97+0- 90 

U - y ' -0.57 

n 99+0 0 3 
u.yy_ 0 04 

1 nn+ 0 - 06 
l.UU-o 05 

1 oi +°- 06 
l-Ul —0.05 

Wan 

Rq 

1 00 +U-iU 
i.uu_ 0 10 

i-oo±S;8| 

0 QQ+9-U8 

U.yy_ 0 08 

0 99+ 0 - 04 

1 Q3+0-07 
l.UO_o 07 

1 00 +O ' 03 
l.UU_o 03 

1 QQ+0-07 
l.UO-o 07 

i-03±g;“ 

Vc(R o) 

apal 5 

A/200 

P200 

C 

1 - 08 -o:i3 

1 30+ 0 - 84 
l.OU —0 49 

1 4^+0.92 
—0.70 

1 1 9+0.20 
± ' ±z -0.23 

1 09+I.I9 
1 - OZ -0.64 

0 98 +(J ‘ Ufci 
u.ho _ 0 04 

0 93+ 0 - 25 
u - yo -0.17 

1 nc+0.51 
4 -UO_ 0 41 

i-02tS:ll 

0.88l°;“ 

0 98 +(J ' U4 
o.yo _ 0 03 

0 90 +0 ' 18 
o.yu _ 0 13 

0 92 4 ” 0 - 32 
U ' yz -0.26 

0 97 + 0 - 10 
- 0.10 

0 89 +0 ' 20 
u.oy _ 0 is 

0-9811- 
0 90 +0 ' 14 

u.yu-o 11 

0 95 +0 ' 30 
u.yo_o 25 

0-98i°;“ 

0-84l°l® 
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M H aio[10 12 Mo] r Ha i 0 [kpc] 

Fig. 6. — Posterior parameter distributions of the three halo 
parameters (scale mass, Mp a i 0 \ scale radius, rpalo ; flattening 
along z-axis, q z \ see Eq. 7). The blue-shaded histograms show 
the marginalized posterior distributions for each parameter, the 
red lines indicate the median values, the dashed lines give the 68% 
confidence intervals, and the black vertical lines mark the true sim¬ 
ulation values. The contour plots show the two-dimensional pos¬ 
terior distributions for parameter pairs. Uncorrelated pairs should 
show roundish contours, whereas correlated parameters like rpalo 
and Mtfaio have elongated distributions. All model parameters are 
recovered within the respective uncertainties. 


for the three respective parameters. The median of 
this distribution is indicated by a red vertical line, and 
the 68 % confidence intervals are marked by dashed 
vertical lines. The true parameter value from the analog 
simulation is shown by a black solid line. Ideally, the 
red line should lie on top of the black. 

According to Fig. 4 we can recover an underlying 
gravitational potential with a Pal 5-like stream to an 
accuracy of ~ 10% and with an uncertainty of 20-30% 
in each halo parameter when we use the weighted 
overdensity analysis approach. This also holds for 
quantities derived from our model parameters, keeping 
in mind that the parametrizations of the Galactic 
potential in the 7V-body simulation and in the streak¬ 
line modeling is identical. In Fig. 4, we also show 
the (more frequently used) virial mass, M 2005 virial 
radius, r 2 oo> and concentration, c, derived from the 
scale masses and scale lengths of our trial models. 
As expected, the uncertainties on these quantities are 
comparable to the ones of the original model parameters. 

Interestingly, the acceleration at the (assumed) lo¬ 
cation of Pal 5, ap a i 5 , shows a factor of 2 lower spread 
than, e.g., M 200 (Tab. 5). This is due to the fact that 
scale mass and scale radius of the halo are strongly 
degenerate (see lower left panel of Fig. 6 ), and a wide 
range of combinations of these two quantities can 
produce very similar accelerations at the location of 
Pal 5. Reporting accelerations at given points within 
the Milky Way halo, instead of halo masses within a 


certain radius, or virial masses, may therefore be a 
more reasonable way of comparing results from different 
studies - especially if the halo is found to be aspherical, 
since then the term ’enclosed mass’ is ill-defined. 

3.1.2. Palomar5 parameter results (analog simulation) 

For our streakline models, five out of 10 model param¬ 
eters determine the properties of each trial cluster. The 
heliocentric distance of Pal 5, dp a i 5 , and its two proper 
motion components, /i a cos (5) and pis, set its orbit. The 
present-day mass, Mp a i 5 , and mass loss rate, dM/dt , fix 
its mass evolution. All of these cluster parameters are 
recovered within the 68 % confidence intervals from all 
four methods (Tab. 5) - the only exception is again the 
interpolated centerline method, which fails at recovering 
the cluster mass. 

From Fig. 4 it is obvious that the three orbital pa¬ 
rameters are significantly better constrained through the 
data than the mass evolution parameters. The weaker 
performance on these mass parameters is most certainly 
due to the weak dependence of the tidal radius (alas 
the key quantity that links the tails and the cluster) on 
the cluster mass (see Eq. 2). Nevertheless, all analysis 
approaches except for the interpolated centerline method 
recover the cluster mass, even though the departures 
can be as large as +123% and —83% if we use just the 
overdensities. Yet, with a precision of about ±50%, 
the weighted overdensities analysis approach provides a 
robust measure for a quantity that is hard to determine 
observationally (see Sec. 4). 

The mass loss rate is only weakly constrained and 
shows an asymmetric posterior distribution with a tail 
towards high mass-loss rates (Fig. 7). Even so, it is 
intriguing that the shape of a tidal stream can give 
insights on the mass evolution of its progenitor. 

The full strength of tidal streams as high-precision 
scales can be seen in the orbital parameters. The 
precision and accuracy of the orbital parameters is 
better than « 10% for all four methods. The weighted 
overdensity analysis approach recovers the distance to 
the mock Pal 5 with an accuracy of 1.00 and a precision 
of ±4%. The proper motion components are also 
accurately recovered by this approach to a precision of 
« 7%. As can be seen in Fig. 7, the two components 
are highly correlated due to relatively close alignment 
of the tails and the cluster orbit, that is, the direction 
of the cluster’s transverse direction of motion is very 
well constrained through the tails, whereas the cluster’s 
transverse speed allows for some variation. 

3.1.3. Solar parameter results (analog simulation) 

Just like for Pal5’s distance and its proper motion, 
the two Solar parameters in our modeling, V ta n and 
Rq, are recovered with remarkable accuracy by all 
four methods. Using the weighted overdensities, the 
precision can be as good as 7% for the Solar trans¬ 
verse velocity, and 3% for the Galactocentric radius of 
the Sun. This clearly demonstrates the power of cold 
tidal streams in constraining length and velocity scales 11 . 

11 As long as they lie within the orbit of the progenitor. 
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Fig. 7.— Same as Fig. 6 but for the cluster parameters (present-day mass of Pal 5, Mp a | 5 ; mass loss rate, dM/dt ; heliocentric distance, 
dpai 5 ; proper motion components, cos(e)) and ^ 5 ). All model parameters from the analog simulation (black lines) are recovered within 
the 68 % uncertainty intervals (dashed lines). 


1 

M 



R O [kpc] Vtan [km/s] 

Fig. 8 .— Same as Fig. 6 but for the Solar parameters (So¬ 
lar transverse velocity with respect to the Galactic center, Vtan ; 
Galactocentric radius, Rq). The Solar parameters from the ana¬ 
log simulation (black lines) are recovered with high accuracy. 


As a sanity check, we computed the circular veloc¬ 
ities at the assumed Galactocentric radius of the Sun, 
Vc{Rq ), for all parameter sets in our MCMC chains. 
All methods recover the expected value with high 


accuracy and precision (Tab. 5). This is reassuring, 
but should be taken with a grain of salt as in both 
7V-body simulation and streakline modeling we used the 
same (fixed) potentials for the disk and the bulge (see 
Sec. 2.4). Moreover, in both cases we used the same 
parametrization for the dark halo potential, which for 
the real Milky Way we unfortunately do not know. This 
should be kept in mind when interpreting the results 
from the observations of Pal 5. 

3.2. Palomar 5 observations 

The posterior parameter distributions of our modeling 
of the Pal 5 observations look basically the same as the 
ones from the analog simulation (see Fig. 10-12), just 
with slightly different median and uncertainty values 
(Tab. 6). This was to be expected since we chose the 
parameters of our analog N -body simulation based on 
preliminary modeling of the observations. Hence, since 
the analog simulation probes a configuration that is 
very similar to the real Pal 5, the underlying systematic 
uncertainties should also be comparable. 

For the weighted-overdensity approach we show 
the median model in Fig. 9. The epicylcic loops of the 
stream stars are clearly visible and the closest loops 
can be associated with the overdensities we found in 
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TABLE 6 

Results from fitting to the observations 



Overdensities 

Overdensities 

Weighted 

Interpolated 

Parameter 


+ radial velocities overdensities 

centerline 

M Halo [10 12 M©] 

1 60 +u ' 72 
l.OU —o 63 

i 7 p:+0.70 

4 - 1 °-0.66 

1 e;q++^ 

l.OO_0 37 

1 +^-4:5 

i '° i -0.43 

rHalo [kpc] 

oq 9 + 21.5 
ooz -13.6 

41-8111.0 

37.91®' g 

37-Olg-g 

Q.z 

1 ni+0.23 

1-U1 —0.24 

o.84l° 0 ;?S 

0.95l°+ 

l-OOlgili 

Mpal 5 [10 3 M©] 

16.81 + 2 

1 8 0 + 12.9 

lo-O_i 0- 6 

ie.olS;g 

14-5111 

dM/dt [10 3 M 0 Gyr" 1 ] 

4-6 ±t:l 

10. ill;! 

7.9l®;3 

1 Q+ 2 ' 4 

4 - y -i.o 

dpal 5 [kpc] 

23.52l}+ 

23.50li;iS 

23.581q'®2 

oq -| o+0. 87 
zo.rz _ 1 00 

lla cos(e)) [masyr -1 ] 

—2 49 + 0 - 25 
Z -^ y -0.28 

—2 40 +0 - 21 
Z>4:U -0.23 

—2 J9 + 0 - 15 
z.oa _ 0 17 

— 2 50 - * -0 ’ 17 
z - ou -0.18 

1*6 [masyr -1 ] 

—2 43+ 0 - 22 
Z -^-0.24 

o qo+ 0.20 

2 .dO 0 21 

-2.36l“; 1 4 

_2 44+ 0 - 15 
z '^-0.16 

Vtan [kms -1 ] 

252.9 

254.5^21.4 

254.3l,g'[ 

251 l+ 18 - 4 

z 0l.l_i7 o 

Rq [kpc] 

Q qr) + 0.37 

o-5^ 0 .42 

o qc+0.33 
O.OO 0 35 

8-301 °oii 

8 29 +0 ‘ 30 
°- zy -o.3i 

V c ( J R 0 )[kms- 1 ] 

241 7+ 4i - 8 
z ^ 4 - ‘ -24.6 

230.6^i";| 

233.0llg;5 

238.2lg 1 3 5 

ipal 5 [pcMyr -2 ] 

2 91+ 1-65 
-0.99 

2 60 +0 ‘ 87 
Z.OU_o 63 

2 «i+ 0 - 56 
Z - Di -0.44 

9 qq+0.53 
^•o9_0.48 

M 200 [10 12 M 0 ] 

4 - ' U -0.84 

i-5»±8:lg 

1 '5fi+ 0 - 42 
l-OO_ 0 42 

1 op:+0.49 
l-OO_ 0 43 

^200 [kpc] 

200.5138-2 
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Fig. 9. — Streakline model (green semi-transparent points) of Pal 5 using the best-fit parameters from the weighted-overdensity analysis 
approach. The present-day position of the cluster is marked by a red star, whereas the locations and uncertainties of overdensities (upper 
panel) and radial velocity stars (lower panel) are shown by black markers. Just like for the Pal 5-analog simulation, the streakline model 
traces the periodic overdensities very accurately close to the cluster. 
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M H aio[10 12 Mo] r Ha i 0 [kpc] 

Fig. 10.— Posterior parameter distributions of the three halo pa¬ 
rameters (scale mass, Mn a lo\ scale radius, rpalo\ flattening along 
z-axis, q z \ see Eq. 7) as recovered from the Pal 5 observations using 
the weighted overdensity analysis approach. The green-shaded his¬ 
tograms show the marginalized posterior distributions for each pa¬ 
rameter, the red lines indicate the median values, and the dashed 
lines give the 68% confidence intervals. Our best-fit NFW halo 
is slightly oblate, has a large scale radius, and a mass of about 
(1.6 ±0.4) x 10 12 M q . 

the SDSS data. At about 5 deg distance from the 
cluster, the epicyclic loops overlap so much that no clear 
substructure is generated, especially since our streakline 
models have no scatter in their escape conditions, and 
hence in a more realistic simulation the structures would 
wash out more easily. 

We will discuss the results from the different anal¬ 
ysis approaches in detail for the three parameter groups 
in the following. 

3.2.1. Halo parameter results (observations) 

Just like for the analog simulation, we get similar 
results from all four analysis approaches, and similar 
to the analog simulation results, the results from the 
observations also get more precise the more constraints 
we use (Tab. 6). Moreover, all results for the halo pa¬ 
rameters from the different approaches agree with each 
other within their respective 68% confidence intervals. 
In the following we will discuss the results from the 
approach using the weighted over densities as it shows, 
on average, the smallest uncertainties. The results from 
this analysis approach are very similar to the results 
using the interpolated centerline. 

Assuming that the dark halo potential of the Milky 
Way has an NFW shape, we get from our modeling of 
the Pal 5 observations that the halo is best represented 
by a scale mass of Mpaio — (1.6 ± 0.4) x 1O 12 M 0 , 
a scale radius of rpaio — 38^g kpc, and a flattening 
along the z-axis of 0.951 q i 2 (Fig- 10)- The respec¬ 
tive virial mass, virial radius and concentration are 
M 2 oo = (1.6 =b 0.4) x 10 12 M o , r 2 oo = 195tjg kpc, and 


cm 5.1 1*1 i' 2 . This is higher than most recent estimates of 
the Milky Way’s virial mass, and also the concentration 
value is significantly lower than most studies find (see 
Sec. 4.1). This may be due to the fact that the scale 
radius of the assumed NFW halo lies well outside Pal5’s 
orbit (apogalactic radius of ~ 19 kpc) and is therefore a 
vast extrapolation out to about 200 kpc, and/or that the 
NFW potential is a poor representation of the Galactic 
potential in the inner 19 kpc, and/or that we fixed the 
disk and bulge potentials. 

A more sensible and also better constrained quan¬ 
tity is the acceleration at the present-day position of 
Pal 5. We find that ap a is = 2 . 61 q '.4 P c Myr -2 , which 

corresponds to (0.8llo'iI) x 10 - 10 ms -2 (Tab. 6 ). Using 
Eq. 9, this corresponds to a local circular velocity of 
Vc (rpais) = 217^9 km s 1 at the Galactocentric radius 
of Pal5 (rpai 5 = 18.81o!t kpc). From the acceleration 
we can also calculate an “equivalent mass”, that is, the 
enclosed mass assuming the mass inside the present-day 
radius of Pal 5 was spherically distributed: 

apai 5 r 2 pal 5 G- 1 = (2.06±37) * 10 U M Q . (16) 

In the symmetry plane of our Galactic potential we 
do not have these problems of quoting sensible values 
for comparison with other studies. A common way to 
express the acceleration or enclosed mass within the 
galactic disk is the circular velocity. Even though we 
seem to limit ourselves by choosing a fixed parametriza- 
tion for the dark halo, we find that all analysis 
approaches produce median posterior values for the 
circular velocity at the Solar radius, Vc{Rq ), between 
231 km s -1 and 242 km s -1 . That is, all recovered halo 
parameters produce reasonable halos that are in agree¬ 
ment with most estimates for the circular velocity at the 
Solar circle (see Sec. 4 for a more detailed comparison to 
other observational constraints). This is reassuring and 
tells us that our choice of parametrization is at least not 
obviously wrong. The weighted overdensity approach 
yields a circular velocity of Vc(Rq) = 233.0t}o!o kms -1 . 

As we will see in the following, our results for the 
other model parameters also agree with independent 
estimates from the literature. 

3.2.2. Palomar 5 parameter results (observations) 

The posterior parameter distributions for the 5 
cluster parameters (Fig. 11) look remarkably similar in 
extent and shape to the posterior distributions from 
the analog simulation (cf. Fig. 7). With the weighted 
over density analysis approach, we infer a cluster mass of 
M Pal5 = (16.0±|;|) x 1O 3 M 0 . 

As with nearly all other parameters, the four ap¬ 
proaches produce similar results, which all agree within 
their 68% confidence intervals with each other (Tab. 6). 
One exception is the mass loss rate, for which the 
interpolated centerline method yields a comparatively 
low value with small uncertainties. We haven’t observed 
such a bias in this particular parameter for the analog 
simulation, hence we cannot relate this discrepancy to 
any known effect. Looking at the posterior distribution 
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|i 6 [mas/yr] Mp a i5[10 3 M o ] dM/dOMoyr 1 ] d Pa i 5 [kpc] |i a cos(5) [mas/yr] 


Fig. 11.— Same as Fig. 10 but for the cluster parameters (present-day mass of Pal 5, Mp a i 5 ; mass loss rate, dM/dt\ heliocentric distance, 
dpai 5; proper motion components, /j,a.cos(5) and ^5). 


TABLE 7 

Orbital parameters of the best-fit model of each analysis approach 



Overdensities 

Overdensities + radial velocities 

Weighted overdensities 

Interpolated centerline 

Present-day distance 

18.70 kpc 

18.68 kpc 

18.77 kpc 

18.34 kpc 

Apogalactic distance 

19.06 kpc 

19.07 kpc 

19.15 kpc 

18.67 kpc 

Perigalactic distance 

8.30 kpc 

7.83 kpc 

7.54 kpc 

7.97 kpc 

Orbital eccentricity 

0.39 

0.42 

0.43 

0.40 


of this generally weakly constrained parameter, however, 
we find that it tends to be asymmetric and bimodal 
(Fig. 11). We already pointed this out for the Pal 5 
analog, but here the bimodality appears to be more 
pronounced. The interpolated centerline approach seems 
to strongly prefer the peak with the lower mass loss 
rate. The other analysis approaches yield larger values 
for dM/dt but also larger uncertainties. 

The proper motion in right ascension shows a slightly 
broadened distribution with respect to the analog sim¬ 
ulation, even though the same parameters also showed 
a lower kurtosis for the mock Pal 5. For the observed 
Pal 5 the posterior distribution in fi a cos(£) may even be 
bimodal. Since cos (5) and are strongly correlated, 
the proper motion in declination also shows a broader 
distribution than we got for the analog simulation. The 
median proper motion values and the corresponding 68% 


confidence intervals are /i a cos(5) = — 2.39^q!it mas Y r 1 
and ns = — 2.36^q’.i 5 m asyr -1 , which corresponds to a 
precision of better than 7%. 


The other two cluster parameters are well behaved. 
Pal5’s heliocentric distance we recover with a remark¬ 
able precision of 3% to be 23. 6 ^ 0 !? kpc. As for the mock 
Pal 5, we find that the geometry of Pal 5’s tidal stream is 
very powerful in constraining length and velocity scales. 


It is reassuring that the four different analysis ap¬ 
proaches for comparing the streakline models to the 
observational data points yield comparable results. 
The models using the median values of the posterior 
parameter distributions from the four approaches are 
hardly distinguishable from each other, and also the 
orbital parameters of the median models all fall within 
a narrow range (Tab. 7). 
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R© [kpc] Vtan [km/s] 


Fig. 12. — Same as Fig. 10 but for the Solar parameters (So¬ 
lar transverse velocity with respect to the Galactic center, Vtan ; 
Galactocentric radius, Rq). 

3.2.3. Solar parameter results (observations) 

As expected from the modeling of the analog simu¬ 
lation, the two parameters of our streakline models that 
describe the location and velocity of the Sun with respect 
to the Galactic rest frame, and hence to Pal 5, turn out 
to be very well constrained. For the transverse velocity 
of the Sun with respect to the Galactic center we get 
Vtan = 254 d= 16kms -1 , and for the Galactocentric dis¬ 
tance of the Sun we get Rq = 8.301q!25 kpc. Again, the 
four analysis approaches yield statistically indistinguish¬ 
able results. Moreover, these values (and others of the 10 
model parameters) agree very well with independent de¬ 
terminations from the literature. These will be discussed 
in the following. 

3.2.4. Hyper-parameter results (observations) 

In Sec. 2.6 we described the four hyper-parameters of 
our modeling process, that is, the integration time, the 
smoothing in positional uncertainty of the over densities, 
the smoothing in the positional uncertainties of the radial 
velocity stars, and the smoothing in the radial velocities. 
These were allowed to vary between zero and nearly arbi¬ 
trarily large values. However, the posterior distributions 
for these hyper-parameters turn out to be absolutely rea¬ 
sonable and well-behaved. The distribution of integra¬ 
tion times peaks at 3.4+oJ Gyr, telling us that the parts 
of the Pal 5 stream for which we detect significant over¬ 
densities within the SDSS footprint is of about this age. 
The positional uncertainty of the overdensities peaks at 
(0.13 zb 0.03) deg, whereas the positional uncertainty of 
the radial velocity stars lies at 0.21 q.2 deg. Finally, the 
smoothing in radial velocity is (0.7 zb 0.3) kms -1 . 

4. DISCUSSION 

As we have shown in the previous Section, the four 
different approaches for comparing the streakline models 
to the observational data yield results that are generally 
consistent among each other. This gives us confidence in 
the stability of the results, that is, the individual pos¬ 
terior parameter distributions for each model parameter 
are not artifacts of a single fitting method. To ensure 
that the analysis approaches do not fail collectively, we 
generated a simulation of a Pal 5-analog on a similar orbit 



TABLE 8 

Mass of the 

Milky Way based on Pal 5 

Measure 

“Weighted overdensities” estimate 

M 200 

(1.69 ±0.42) x 10 12 Mq 

M(r < 100 kpc) 

(0.90 ± 0.20) x 10 12 M 0 

M(r < 50 kpc) 

(0.43 ±0.11) x 1O 12 M 0 

M{r < r p a i 5 ) disk 

2 - 14 to;|| x 10 11 M© 

UPal 5 

2.61±°;®s pcMyr -2 

Vc(rpal 5) 

217+ 2 g kms -1 

a Pal 5 r p a l 5 G — 1 

2Mt° 0 A 3 7 x 1011 M © 


and in a similar potential as the real Pal 5. Our method¬ 
ology passed this test in that it successfully recovered 
all model parameters. Yet, the analog TV-body simula¬ 
tion is quite idealized compared to the real Pal 5: in the 
simulation we know for example the exact shape of the 
underlying galactic potential, and we fixed the mass and 
extent of the disk and bulge to some reasonable values. 
In order to test if the outcomes of the Pal 5 modeling is 
reasonable and close enough to reality, we have to cross¬ 
check our results with estimates from other, independent 
methods. This is what we will do in the following for the 
three model parameter groups. 

4.1. Comparison of halo parameters to literature values 

Weighing the Milky Way and its components has been 
one of the major goals in Galactic dynamics within the 
last few decades. Since the total mass of the Galaxy 
is vastly dominated by the extended dark halo, weighing 
the entire Galaxy basically boils down to determining the 
mass profile or gravitational potential of this halo. Var¬ 
ious attempts have been made at constraining its virial 
mass, or the halo mass enclosed within a certain limiting 
radius, using all kinds of available tracers of the Galac¬ 
tic matter distribution. Some approaches also attempted 
to constrain the shape of the dark matter distribution, 
i.e. its flattening or triaxiality. Here we list a few of the 
more recent measurements, first for the enclosed mass 
and then for the shape of the dark halo. 

4.1.1. Enclosed mass 

Most of the estimates listed in the following are 
summarized as a rotation curve in Fig. 13, in which we 
collected enclosed mass estimates and circular velocity 
estimates from these references. However, we note that a 
comparison between the various studies is difficult due to 
the different measures quoted by the authors. Moreover, 
different tracers probe different radii, and some probe 
the Galactic potential in the MW disk, whereas others 
probe it at different points within the (not necessarily 
symmetric) halo. In an attempt to facilitate comparison, 
we compiled a list of “classical” measures from our own 
modeling of Pal 5 (Tab. 8), but would like to point out 
that Pal 5 essentially probes the Galactic mass within 
its apogalactic radius of ~ 19 kpc. Our estimate for the 
MW’s mass beyond that radius or its virial mass solely 
relies on the assumption that the dark halo potential is 
of NFW form. From our modeling there is no indication 
that this is the correct parametrization for the dark 
halo potential. In fact, most of the most recent studies 
indicate that the virial mass of the Milky Way halo is 
lower than our best-fit value. This may indicate that 
the NFW potential parametrization is a poor choice in 
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Fig. 13.— Collection of circular velocity estimates from the literature. Our estimate for R « 19kpc bridges the gap between disk 
tracer estimates and halo tracer estimates. The black solid lines give the 68% confidence intervals for the rotation curve resulting from our 
modeling of Pal 5. The NFW parametrization that we chose for the dark halo potential does not seem to fall off quickly enough at radii 
larger than 19kpc to reproduce results from other investigations. References: W99 (Wilkinson & Evans 1999); S03 (Sakamoto et al. 2003); 
B05 (Battaglia et al. 2005); S07 (Smith et al. 2007); X08 (Xue et al. 2008); G10 (Gnedin et al. 2010); K10 (Koposov et al. 2010); W10 
(Watkins et al. 2010); Mil (McMillan 2011); B12 (Bovy et al. 2012); D12 (Deason et al. 2012); K12 (Kafle et al. 2012); B14 (Bhattacharjee 
et al. 2014); G14 (Gibbons et al. 2014). 


the case of the Milky Way. 

Irrgang et al. (2013) fit three different, simple an¬ 
alytic galaxy models, each consisting of a disk, a 
bulge and a (spherical) halo, to a range of dynamical 
constraints. They compute various quantities for 
their models, of which we mention here only the mass 
enclosed inside a radius of 50 kpc, as it involves the 
least extrapolation. Depending on the parametrization 
they choose, they get values for M{r < 50 kpc) between 
(0.46 ± 0.03) x 10 12 M e and (0.81^;^) x 1O 12 M 0 . 
It is noteworthy that these estimates exclude each 
other, and solely differ in the choice of parametrization. 
We therefore focus on our mass estimate for Pal5’s 
apogalactic radius of about 19 kpc, and only give the 
other mass estimates for completeness. 

However, since Pal 5 is located well above the disk, 
it is outside any symmetry plane, which makes a 
definition of an enclosed mass difficult. A better 
estimate can be made, if we instead calculate the 
acceleration vector on Pal 5 at its current position from 
the analytic potential. For its absolute value, we get 
apai 5 = 2 . 61^^44 pcMyr -2 ( 0 . 8 lt°;i 4 x 10 -lo ms -2 ), 
which we can convert into a “circular velocity” 
(using Eq. 9) of V c (r Pa i 5 ) = 217± 2 gkms _1 , 

or an equivalent spherically distributed mass 
yielding ap a Z 5 7 p ai5 G _1 = 2.06±° 0 ?£ x 1 O U M 0 . 

If we assume that Pal 5 lies in the plane of 
the Galactic disk at this radius, we arrive at 
M(r < r P aih)disk = 2. 141^34 x lO n M 0 . As can 


be seen in Fig. 13, our circular velocity estimate bridges 
a gap between estimates from halo tracers and estimates 
from disk tracers (Fig. 13). It fits well into the general 
trend of the Milky Way rotation curve when combined 
with the estimates listed below. 

The brightest tracers in the Milky Way halo, and 
hence the most obvious targets for Milky-Way-mass 
studies, are the Galactic satellites, i.e. globular clusters 
and dwarf galaxies, but also bright stellar tracers for 
which reasonable distances can be estimated, e.g. blue 
horizontal branch stars (BHB). Results from such 
studies differ widely due to ambiguities in the assumed 
anisotropy of the tracers. Wilkinson & Evans (1999) 
used positions, distances, radial velocities and a few 
existing proper motions of 27 satellites beyond 20 kpc 
to estimate a Milky Way mass of 5 . 4^3 g x 10 11 M 0 
within a radius of 50 kpc. Similarly, Watkins et al. 
(2010) get masses for the Milky Way inside a radius of 
100kpc of M(r < 100kpc) = (0.2 — 2.3) x 1O 12 M 0 , 
depending on the assumed orbital anisotropy of the 
satellites. Their estimate using the anisotropy provided 
by the (limited) data itself is (1.4 db 0.5) x 1O 12 M 0 . 
Sakamoto et al. (2003) use satellite galaxies, globular 
clusters and horizontal branch stars to constrain the 
mass within the inner 50 kpc of the Galaxy to be 
M(r < 50 kpc) = 5 . 3 lo '4 x 1O 11 M 0 . Battaglia et al. 
(2005) use similar data to derive an enclosed mass within 
120kpc of M(r < 120kpc) = 5 . 4+ 2- 4 x 1O 11 M 0 . Smith 
et al. (2007) use RAVE data to measure a mass for the 
inner 50 kpc of M(r < 50 kpc) = 4 . 0 to 7 6 x 10 11 M 0 
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assuming an NFW halo. Analyzing a sample of radial 
velocity measurements for a large sample of BHB stars 
in the halo, Xue et al. (2008) derive a mass within 60 
kpc of (4.0 ± 0.7) x 10 11 M 0 . Gnedin et al. (2010) apply 
spherical Jeans models to a sample of several hundred 
BHB and RR Lyrae stars with distances between 
25 and 80 kpc. They determine an enclosed mass of 
M(r < 80 kpc) = 0.69^012 x 1O 12 M 0 . In a similar way, 
with dynamical modeling of >4000 BHB stars in the 
halo with radial velocities from SDSS/SEGUE, but also 
with allowing for non-sphericity of the tracer distribution 
function, Deason et al. (2012) infer a mass within 50 kpc 
of (4.2 =b 0.4) x 10 n M o , whereas Kafle et al. (2012) 
arrive at a Galactic mass within 25 kpc of 2.1 x 10 n M o 
(without quoted uncertainty). Most recently, Gibbons 
et al. (2014) modeled the Sagittarius stream and arrive 
at M(r < 100kpc) = (4.1 d= 0.4) x 10 11 M 0 , which is at 
the lower end of the spectrum of mass estimates. 

In Fig. 13, we also include some mass estimates, 
or circular velocity estimates, for the Solar environment 
to get a more complete picture of the Galactic rotation 
curve. Koposov et al. (2010) use the GD-1 stream to con¬ 
strain the circular velocity at a distance of 8.4 kpc to be 
(229T18) kms -1 . Bovy et al. (2012) use APOGEE radial 
velocity data on disk stars to constrain the local circular 
velocity at the distance of the Sun from the Galactic 
center to be Vc(R = 8.1 kpc) = (218 db 6) kms -1 . 
McMillan (2011) uses a range of tracers and constraints 
to arrive at Vc{R = 8.3kpc) = 239.2 db 4.8kms -1 . Our 
estimate for the circular velocity at the Solar radius is 
Vc(Rq) = 233.01 1 q o kms -1 , an d f° r Solar radius 
we get Rq = 8.30!o‘. 25 kP c (see Sec. 4.3). 

Ultimately, tracers in the disk and in the halo have to 
be combined and modeled altogether. Following such an 
approach, Bhattacharjee et al. (2014) infer a mass within 
200kpc of M(r < 200kpc) = (0.68 ± 0.41) x 1O 12 M 0 . 
They argue that the uncertainties for the Solar galacto- 
centric distance and the circular velocity at the location 
of the Sun are the two main causes of uncertainties for 
mass estimates in the inner MW, whereas the missing 
constraints on the radial anisotropy of the tracers’ 
phase-space distribution in the halo is responsible for 
the uncertainties at large radii. 

In fact, the mass estimates from the halo tracer 
studies differ widely, and it has to be assumed that 
they suffer from the complexity, i.e. non-sphericity, 
of the Galactic mass distribution, and of streaming 
motions among halo stars. From SDSS data, Bell et al. 
(2008) find the stellar halo to be significantly oblate, 
but furthermore show that a significant fraction of the 
halo stars (>40%) are in spatial substructures. This 
observation is evidence that large part of the stellar 
halo is formed from infalling satellites (e.g., Bullock & 
Johnston 2005). It is unclear what kind of biases such 
coherent structures could cause in the modeling of the 
stellar halo. 

4.1.2. Shape of the dark halo 

While the coherent debris of Galactic satellites is 
unpleasant noise for the modeling of the stellar halo, 


it may offer the most powerful way of measuring the 
mass distribution of the Milky Way. From our point 
of view from within the Galaxy, tidal streams like the 
Sagittarius stream (Ibata et al. 1994; Dohm-Palmer 
et al. 2001; Majewski et al. 2003; Belokurov et al. 2006b, 
2014) or the GD-1 stream (Grillmair & Dionatos 2006b) 
span huge arcs on the sky, and as coherent phase-space 
structures simultaneously probe large parts of the Milky 
Way potential without the need for a tracer density 
profile or assumptions on orbital anisotropy. Moreover, 
the hope is that debris streams are also sensitive to the 
shape of the mass distribution and therefore can give 
us insights to the 3D mass profile of the Milky Way’s 
dark halo. However, the many studies on the Sagittarius 
(Sgr) stream have shown that inferring the shape of the 
dark matter halo from a stream as complex as Sgr is 
non-trivial. 

Since the first models of the disrupting Sgr satel¬ 
lite by Johnston et al. (1995) and Velazquez & White 
(1995), the interpretation of its debris has led to many 
different conclusions. Early numerical models of the Sgr 
stream were made in spherical halo potentials, and the 
quantity and quality of data on the discovered debris 
was not conclusive enough to pinpoint a specific shape 
of the halo. Ibata et al. (2001) pointed out that a 
large flattening of the halo potential would lead to a 
significant orbital precession between the leading and 
the trailing tail, and hence they argued against potential 
flattenings with q z < 0.9. When Sgr was traced along 
a nearly perfect great circle on the sky in 2MASS data 
(Majewski et al. 2003) this hypothesis seemed to be 
further confirmed. However, Helmi (2004a) argued that 
the Sgr stream was in fact too young to be sensitive 
the halo shape, and that the halo potential could 
therefore be either oblate, spherical or prolate. With 
more kinematic data at hand, Helmi (2004b) finally 
argued for a prolate halo shape, while Law et al. (2005) 
and Johnston et al. (2005) argued for an oblate halo to 
explain the small degree of orbital precession observed 
in the Sgr debris. With successive detections of new 
parts of Sgr’s tidal debris, this discussion has yet become 
more complicated. Moreover, Belokurov et al. (2006b) 
and Koposov et al. (2012) have found bifurcations in 
Sgr’s leading and trailing stream, which are still not 
understood. Fellhauer et al. (2006) and Penarrubia 
et al. (2010) tried to explain these new features with 
multiple wraps of the stream around the Milky Way, or 
with rotation in the Sgr progenitor, respectively. Yet, to 
the present day there has been no conclusive model that 
can explain all observed Sgr features. 

Using a large set of collisionless TV-body simula¬ 
tions, Law & Majewski (2010) reproduced many parts of 
the observed debris. However, they found the necessary 
shape of the dark matter halo to be rather peculiar: 
the observed positions, distances and radial velocities 
of Sgr debris stars can be explained by a satellite 
orbiting within a triaxial halo potential that has its 
minor axis pointing towards the Sun within the disk. 
That is, the dark matter distribution resembles an 
oblate halo perpendicular to the Galactic disk, which 
we see face on when we look towards the Galactic 
Center. Deg & Widrow (2013) confirmed these findings 
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using also other observational constraints and a dif¬ 
ferent form of parametrization for the Galactic potential. 

However, several problems have been pointed out 
with such a configuration. Debattista et al. (2013) tried 
to produce a stable galaxy configuration resembling the 
proposed scenario of Law & Majewski (2010) using a live 
dark matter halo, and found it to be highly unstable. 
Moreover, Pearson et al. (2015) demonstrated that the 
thin morphology of the Pal 5 stream cannot be repro¬ 
duced in such a triaxial configuration. Such cross checks 
of stream modeling results are obviously vital; a certain 
Galactic potential that can produce one specific stream 
also has to be able to yield models for all other streams. 
Interestingly, Pearson et al. (2015) found that it is easily 
possible to find a stream model that has the long, thin 
and curved morphology of the observed Pal 5 stream in 
a spherical halo potential. Moreover, the authors point 
out that a model, which fits the morphology within 
a spherical halo potential, also naturally reproduces 
Pal5’s radial velocity gradient, making a strong case 
for a nearly spherical dark halo potential. Ibata et al. 
(2013) have in fact shown that it is possible to reproduce 
the Law & Majewski (2010) observational constraints 
in a galaxy potential that is spherical. However, they 
found the necessary potential to be peculiar as well, 
in that it has to have a rising rotation curve at large radii. 

From cosmological (dark-matter-only) simulations 
we expect dark matter density profiles to be triaxial, 
where the expected axis ratios lie between 0.4 and 0.6 
for the short, intermediate, and long axes (Dubinski 
& Carlberg 1991; Allgood et al. 2006). While these 
axis ratios should be prevalent at large galactocentric 
radii, Dubinski (1994) and Debattista et al. (2008) 
showed that condensations of baryons leads to rounder 
halos in centers of galaxies. Vera-Ciro & Helmi (2013) 
therefore suggested a potential form that is triaxial 
in the outer parts and more spherical in the inner 
parts. They showed that with such a configuration 
and an inner potential flattening of about 0.9 it is 
possible to reproduce a Sgr configuration like the Law 
& Majewski (2010) model. A flattening/triaxiality that 
varies with Galactocentric radius may in fact be a way 
to resolve many conflicts between existing Milky Way 
mass estimates. 

Yet, Sgr’s story is not over: new observational data 
from Belokurov et al. (2014) show that Sgr’s tidal debris 
goes in fact out to 100 kpc. This debris structure cannot 
be reproduced by the Law & Majewski (2010) model. 
While the halo configuration of Gibbons et al. (2014) 
(Milky Way halo with a spherical but sharply cut-off 
density profile) can reproduce the apocentric distances 
of observed Sgr debris, and also the precession of Sgr’s 
orbital plane, it has not been shown to reproduce all 
available observational features. Thus, there is still 
no model of the Sgr satellite that can reproduce all 
observational constraints simultaneously. 

The only other stream that has been modeled to 
infer the shape of the Milky Way’s dark halo has used 
the 63 deg long GD-1 stream (Koposov et al. 2010; 
Bowden et al. 2015). The GD-1 stream is very thin but 


lacks any obvious progenitor satellite, hence it is most 
probably the stream of a dissolved globular cluster. 
In contrast to the Sgr satellite it orbits in the inner 
halo. Koposov et al. (2010) estimate the pericenter 
and apocenter of its orbit to be 14 kpc and 26 kpc, 
respectively. By fitting orbits to sky positions of the 
debris, distances along the stream, radial velocities, 
and proper motions, the authors conclude that the 
best-fitting orbits constrain the dark halo to be nearly 
spherical. They place a lower limit on the flattening of 
the dark halo potential, saying that with 90% confidence 
the halo flattening is q z > 0.89. Recently, Bowden et al. 
(2015) confirmed these results with more sophisticated 
streakline modeling of the stream. 

Our results for the flattening also suggest that the 
inner 19 kpc of the Milky Way’s dark halo are rather 
spherical (assuming that our choices for the disk and 
bulge potentials were reasonable). The results from 
the different fitting approaches all agree within their 
uncertainties with each other, and they are actually 
all compatible with a spherical halo. Our best esti¬ 
mate is q z = 0.95lo'i2? which supports the notion of 
Vera-Ciro & Helmi (2013) of having a slightly oblate 
inner halo, and a possibly less relaxed and probably 
triaxial halo at larger radii. On the other hand, our 
result is also compatible with the nearly spherical 
“phantom” halos predicted by alternative gravity the¬ 
ories (Famaey & McGaugh 2012; Liighausen et al. 2014). 

Yet, the existence of the Sgr dwarf galaxy and the 
two Magellanic Clouds should remind us that the halo 
is not necessarily a relaxed structure, but a constantly 
growing and probably highly substructured potential 
well. This is also suggested by Saha et al. (2009) based 
on observations of asymmetries in the Galactic H I disk. 
In fact, the assumptions of steadiness and symmetry 
may be the main reasons for some of the discrepancies 
pointed out above and may pose the philosophical ques¬ 
tion of what we actually measure when we model tracers 
in the Galactic halo. By applying the FAST-FORWARD 
method to 256 globular cluster streams, which formed 
in a cosmologically motivated, and time-evolving dark- 
matter halo of a Milky-Way sized galaxy, Bonaca et al. 
(2014) demonstrated that we can recover the present-day 
mass and shape of the dark matter halo to an accuracy 
of 5-20%. This, however, only holds for an ensemble of 
streams; single streams can show mass biases of up to 
50%. In order to avoid reporting such biased results, 
detailed comparison to simulated streams on similar 
orbits in similar potentials, and comparison to results 
on the same object from other, independent methods 
should be taken into consideration. Our positive result 
from the Pal 5-analog simulation (Sec. 3.1), and the 
agreement of our fitting results for all other model 
parameters with estimates from completely independent 
methods (as shown in the following Sections) give us 
this necessary confidence in our results out to Pal5’s 
apogalactic radius of 19 kpc. 


4.2. Comparison of Pal 5 parameters to literature 
values 


24 


We use five parameters to model Pal 5 and its orbit. 
The parameters are its distance to the Sun, its proper 
motion (2 components), its mass, and its mass loss rate. 

4.2.1. Heliocentric distance 

Of all the cluster parameters in our modeling, the 
distance from the Sun to Pal 5 is the observationally 
best determined parameter. Dotter et al. (2011) use 
deep HST/WFPC2 data to fit isochrones to Pal5’s stel¬ 
lar population, and determine a distance of 23.55 kpc. 
This estimate comes without an uncertainty, but is sup¬ 
posed to be accurate to within Tlkpc (A. Dotter, pri¬ 
vate communication). The agreement with our estimate 
of 23.58^0.72 kpc is remarkable. Harris (1996) lists a dis¬ 
tance of 23.2 kpc, which is based on CCD photometry by 
Smith et al. (1986). This older value is also in agreement 
with our result. 

4.2.2. Proper motion 

The proper motion of Pal 5 is relatively unconstrained. 
Odenkirchen et al. (2002) point out two different 
measurements coming from photographic plates with 
baselines of four decades each. The estimates for 
(/r a cos(£), pis) are (—1.0, — 2.7) masyr -1 (Scholz et al. 
1998), and (—2.55, —1.93) masyr -1 (Schweitzer et al. 
1993, revised by K. Cudworth in 1998). Our best-fit 
values are (—2.39lo'.i7> —2.36^o!i5) masyr -1 , and lie 
well within the range of the two measurements. Future 
proper motion measurements with HST will test this 
parameter more accurately. 

Pawlowski et al. (2012) point out that Pal 5 and 
its stream lie within the vast polar structure, which is 
a coherently rotating disk of satellites and streams that 
contains most of the Galactic satellites. Our proper 
motion estimate verifies the notion of Pawlowski & 
Kroupa (2014) that Pal 5 is in fact counter orbiting 
within this structure. 

4.2.3. Mass 

Pal5’s mass is not well determined. The Harris 
(1996) catalogue lists an absolute V-band magnitude 
of My = —5.17 mag, which corresponds to L/Lq = 
Iq( m v,o-m v )/2.5 _ 10 000. A possible mass-to-light ratio 
between 0.5 and 3 allows a wide range of cluster masses. 
Our estimate of I6.OI5 9 x 10 3 M© is well within this 
range. However, Odenkirchen et al. (2002) arrive at an 
absolute magnitude of only My = —4.77 =b 0.20 mag, 
i.e. a luminosity of 6 9001/©. They also measure the ve¬ 
locity dispersion of Pal 5 and estimate that their sam¬ 
ple is contaminated by a large fraction of binaries (24- 
63%). After correcting for this, they arrive at a cluster 
velocity dispersion between 0.12 and 0.41 kms -1 , corre¬ 
sponding to a low cluster mass of 5 — 8 x 10 3 M©, when 
fit with a dynamical King (1963) model. This analy¬ 
sis relies on the assumption that the contribution from 
the binaries can be clearly separated from the rest of 
the cluster stars. Uncorrected, the velocity dispersion is 
(1.1 dz 0.2) kms -1 , in agreement with a recent determi¬ 
nation by Kuzma et al. (2015) of (1.3 =b 0.4) kms -1 . Us¬ 
ing a simple virial argument and Pal5’s half-light radius 
of rh ~ 20 pc, we can estimate the virial mass of Pal 5 


to be Mp a r |5 = Krhcrl D G 1 = n x 16 900 M©, where 
(J^d = V^cr is the three-dimensional velocity dispersion, 
and k is some proportionality constant, which relates 
the cluster’s half-light radius to the virial radius. For 
a Plummer model n ~ 1.3, which would yield a mass of 
about 22 000 M©. However, this calculation assumes that 
Pal5’s density profile is well approximated by a Plum¬ 
mer sphere, that mass follows light in the cluster, and 
that the cluster is in virial equilibrium. All three of 
the assumptions may not be fulfilled, since the cluster’s 
surface brightness profile may be better approximated 
by a low-concentration King model (Odenkirchen et al. 
2002), the cluster itself is heavily mass segregated and its 
dark remnant content is unknown (Koch et al. 2004), and 
the cluster may have been tidally shocked recently and 
brought out of equilibrium (Dehnen et al. 2004). The 
fact that we provide an independent measure on Pal5’s 
mass, which does not depend on assumptions about the 
stellar mass function, dark remnant content, binary frac¬ 
tion and virial state, may help to understand the other 
mass estimates better. With more kinematic data on the 
cluster stars and more sophisticated dynamical modeling, 
our result could be easily tested. 

4.2.4. Mass loss rate 

The mass-loss rate of Pal 5 has been estimated by 
Odenkirchen et al. (2003) using SDSS star counts and 
simple assumptions on the drift rate of stars along the 
tails. Using their lower present-day mass for Pal 5 of 
6 000M©, the authors estimate a rate of 4.9M©Myr -1 
from the trailing tail, and 3.8M©Myr -1 from the lead¬ 
ing tail star counts. Our weakly constrained estimate 
of 7.9tg; 3 M© Myr -1 is compatible with these results. 
When we use their cluster-mass independent estimate 
of — M / Mp a i 5 = (0.7 dk 0.2) Gyr -1 and our present- 
day mass of 16000M©, we get dM/dt = (11.2 db 

3.2) M© Myr -1 , which is on the upper end of our esti¬ 
mate. Like for the cluster mass, these estimates are com¬ 
pletely complementary, as our result does not depend on 
star counts or assumed drift rates, but is a free parame¬ 
ter in the modeling process. Deeper imaging of the tidal 
tails will help to assess the mass within the stream and 
estimate the mass loss rate. 

4.3. Comparison of Solar parameters to literature 
values 

The two Solar parameters in our modeling are the 
Galactocentric distance of the Sun, 72©, and its tangen¬ 
tial motion with respect to the Galactic Center, V t an- 
Independent estimates for these quantities exist from var¬ 
ious other methods. 

4.3.1. Galactocentric distance 

The distance of the Sun from the Galactic Center has 
long been unknown to about 10% uncertainty. With 
high-precision measurements of, e.g., the parallax of Sgr 
A*, the orbital motion of the central S-stars, the dy¬ 
namics of the nuclear star cluster, or of maser sources 
in the Milky Way disk, our knowledge on this quantity 
has increased to better than 3% uncertainty within re¬ 
cent years (Fig. 14). Genzel et al. (2010) summarize the 
many different direct and indirect estimates for the dis¬ 
tance of the Sun to the Galactic Center prior to 2010 
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Fig. 14.— The distance of the Sun to the Galactic center, Rq, 
as estimated by recent studies. Our estimate is given on the far 
right. It agrees well with the weighted mean of all measurements 
shown as green solid line with green dashed lines indicating the 
weighted sample variance. For comparison, the IAU standard of 
8.5 kpc is shown as a red dotted line (Kerr & Lynden-Bell 1986). 

to be Rq = (8.25 =b 0.19) kpc, where the uncertainty 
is the variance of the weighted mean of all the mea¬ 
surements. McMillan (2011) fit dynamical models to 
various photometric and kinematical data to arrive at 
Rq = (8.29 zb 0.16) kpc. Similarly, just by sophisticated 
modeling of existing observational data, Schonrich (2012) 
gets Rq = (8.27±0.29) kpc. More recently, Dekany et al. 
(2013) use distance data on thousands of RR Lyrae stars 
in the bulge to get Rq = (8.33±0.05±0.14) kpc. Do et al. 
(2013) use anisotropic, spherical Jeans models to model 
3D kinematic data of the nuclear star cluster, and arrive 
at Rq = (8.92 ±0.58) kpc, and Rq = 8.46lo‘.38 kpc when 
they include data on the orbit of the star SO-2 around 
Sgr A*. Chatzopoulos et al. (2015) improve upon the 
Do et al. (2013) results by using in addition 2 500 line- 
of-sight velocities, 10 000 proper motions, and 200 maser 
velocities. Moreover they also fit the surface density pro¬ 
file of the nuclear star cluster and allow for rotation. 
Their modeling yields Rq = (8.30 ± 0.09 ± 0.1) kpc, and 
(8.36 ± 0.11) kpc if they include the S-star orbits from 
Gillessen et al. (2009). Finally, Reid et al. (2014) de¬ 
termine Rq = 8.34 ± 0.16 kpc from trigonometric par¬ 
allaxes of star forming regions in the Galactic disk. All 
these recent estimates agree very well with our result of 
Rq = 8.30^0.25 kpc. Moreover, our estimate is in good 
agreement with the weighted mean of all these measure¬ 
ments listed above, which is (8.32 ± 0.06) kpc (Fig. 14). 

4.3.2. Solar transverse velocity 

Reid & Brunthaler (2004) determined the tangential 
motion of the Sun with respect to the Galactic Center 
from VLBI measurements of the position of Sgr A* to 
be 0 tan = 6.379 ± 0.026 mas yr -1 in the Galactic plane, 
corresponding to 30.24 ± 0.12kms -1 kpc -1 . Our esti¬ 
mate for this quantity is 0 ta n = 6.43 ± 0.46 mas yr -1 
(30.6 ± 2.2 km s -1 kpc -1 ) and hence in good agreement 
with this strong constraint. Most estimates of the So¬ 
lar tangential motion, Vtam use this measurement of the 
proper motion of Sgr A* and differ only in the way the 
distance to the Galactic Center is estimated. Schonrich 


(2012) uses local kinematics of stars from SDSS DR-8 in 
combination with the Sgr A* measurement to arrive at 
Vtan = 250 ± 9kms -1 . Bovy et al. (2012) use APOGEE 
spectra to model the local Galactic rotation curve. They 
estimate V tan = 242+3° kms -1 . Most recently, Reid 
et al. (2014) use masers in the Galactic disk to measure 
Vtan = 252.2 ± 4.8 km s -1 independently of the motion 
of Sgr A*. They also derive a proper motion of Sgr A* 
from their data, giving 30.57 ± 0.43 km s -1 kpc -1 . All of 
these measurements are in agreement with our best-fit 
value of V tan = 254 ± 16kms -1 . In future modeling of 
Pal 5 we can use, e.g., the proper motion of Sgr A* as a 
strong prior on 6^ an , and hence get better constraints on 
Rq. For now, it is reassuring to see that we get consis¬ 
tent results with these studies that involve very different 
data. 

5. SUMMARY AND CONCLUSIONS 

We have shown that with sophisticated streakline 
modeling of the tidal tails emanating from the Milky 
Way globular cluster Palomar5, and with very little and 
publicly-available observational data (SDSS & radial 
velocities from the literature), we can tightly constrain 
10 different model parameters characterizing the Milky 
Way’s dark matter halo, the cluster Pal 5 itself, and the 
Sun’s position and velocity within the Galaxy. 

To reproduce the characteristic geometry of the 
Pal 5 stream on the sky, we used the streakline method 
outlined in Kiipper et al. (2012) and Bonaca et al. 
(2014), which generates restricted three-body models of 
globular cluster tidal tails in a computationally efficient 
way and therefore allows for scans of a large parameter 
space. We compared the streakline models to locally 
overdense regions of color-selected Pal 5 stars from 
SDSS DR9 (Ahn et al. 2012), which we detected with 
a Difference-of-Gaussians process. The overdensities 
cover a length of 23 deg on the sky and show a typical 
width of 0.1-0.7 deg. Moreover, we used literature values 
for the radial velocity of Pal 5 (Odenkirchen et al. 2002), 
and of 17 radial velocities of giants lying in projection 
and velocity close to the Pal 5 stream (Odenkirchen et al. 
2009). Using a Markov chain Monte Carlo approach, 
and testing the methodology on a Pal 5-analog TV-body 
stream, we constrain the 10 free model parameters of 
the mock Pal 5 with high accuracies and to precisions 
between 3% and 20%. Similar precisions are achieved 
for the observational data, which are all individually in 
agreement with independent estimates from completely 
different and complementary methods (see Sec. 4), 
confirming the validity of our approach. 

Our main findings can be summarized as follows: 

1. Since Pal 5 orbits within an apogalactic radius of 
19 kpc, our best estimate for the mass of the Milky 
Way is limited to this radius. We estimate the 
enclosed mass to be (2.1 ± 0.4) x 10 n M©, corre¬ 
sponding to an acceleration of ap a i 5 = O. 8 II 0.14 x 
10 -10 ms -2 at the present-day position of Pal 5. 
Moreover, we find the halo to have a flattening 
along the z-axis of 0.951 q 12 (Fig. 10). 

2. Our models constrain Pal 5 to have a present-day 
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mass of Mp a i 5 = I6.OI5J x 10 3 M©, and that it 
has been losing mass in the past few Gyr with 
an average mass-loss rate of dM/dt = I.Q^a x 
10 3 M© Gyr -1 . Its heliocentric distance is dp a i$ = 
23.6^?kpc, anc [ ^he proper motion components 
we constrain to be /i a cos(J) = —2.391 q!i 7 masyr -1 
and /is = —2.36^0 15 mas y r_1 (Fig. 11). 

3. The transverse velocity of the Sun with respect to 
the Galactic center is V ta n = 254 db 16kms _1 , and 
for the Galactocentric distance of the Sun we get 
Rq = 8.30io.25kpc (Fig. 12). 

From the remarkable precision of these measurements 
(given the small amount of data), we conclude that long, 
thin and curved tidal streams like the Pal 5 stream are 
exceptionally sensitive to length and velocity scales. The 
sensitivity to mass and acceleration scales is smaller, 
but given the high uncertainties of other methods on, 
e.g., the mass of the Milky Way, it still ranks among 
the most accurate approaches to mass measurements. 
More observational data, such as an extended imaging 
coverage of the Pal 5 stream beyond the SDSS footprint 
or proper motions from, e.g., Gaia, as well as the usage 
of additional prior information on, e.g., the transverse 
velocity of the Sun from proper motion measurements 
of Sgr A*, will significantly improve the precision and 
accuracy that can be achieved. Combining several 
streams will ultimately unleash the full power of tidal 
streams in this respect (Deg & Widrow 2014). 

By testing different analysis approaches for the 
calculation of the model likelihoods, we find that 
radial velocities along the stream are very helpful for 
constraining model parameters. By adding just 17 
velocities, the gain in precision can be as high as 50% 
for individual parameters (see Tab. 5). We also note 
that deeper imaging of the Pal 5 stream would be 
very helpful, as it would increase the accuracy with 
which locally overdense regions of the stream can be 
discovered, and foreground/background fluctuations can 
be removed. The most pronounced overdensities we 
currently find lie in the trailing tail at a distance of 
about 3 deg from the cluster, and have a significance of 
>10cr (see Tab. 1), but many of the farther detections 
are only slightly above the detection limit. Deeper 
imaging may also help to better characterize the nature 
of these overdensities, and conclusively differentiate 
between epicyclic overdensities and, e.g., overdensities 
produced by past encounters of the Pal 5 stream with 
dark matter subhalos (Ibata et al. 2002; Siegal-Gaskins 
& Valluri 2008; Carlberg 2009; Yoon et al. 2011; Carl- 
berg 2012; Carlberg et al. 2012; Erkal & Belokurov 2014). 
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We have shown that globular cluster streams are 
simple and well-understood stellar structures, and 
that we can model them very accurately. We also 
demonstrated that globular cluster streams show density 
variations, which, in parts, can be ascribed to epicyclic 
motion of tail stars. Thus, cold stellar streams in the 
Milky Way halo are promising objects, not just for 
inferring length, velocity and mass scales within the 
Milky Way, but also for studying the dumpiness of 
the dark halo, and hence for drawing conclusions on 
cosmology solely based on simple stellar dynamics. 
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